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APPLICATION OF MULTIVARIABLE SEARCH TECHNIQUES 
TO STRUCTURAL DESIGN OPTIMIZATION 
By R. T. Jones and D. S. Hague 


SUMMARY 

Multivariable search techniques are applied to a 
particular structural design problem, that of determining 
the minimum weight design for a stiffened cylindrical 
shell subject to multiple load conditions. The cylinder 
is stiffened in both longitudinal and circumferential 
directions. Stability analyses are limited to linear 
bifurcation buckling theory. A variety of multivariable 
search techniques embodied in an existing non-linear 
optimization code, AESOP, are applied to this design 
problem. These techniques include elementary single 
parameter perturbation methods, organized search such 
as steepest-descent , quadratic, and FI etc her-Powel 1 
methods, randomized procedures, and a generalized search 
acceleration technique. Design variables are seven in 
number and define stiffener spacings, stiffener dimen- 
sions, and skin thickness. The relative efficiency of 
the techniques are compared. It is shown that a combi- 
nation of search strategies may be superior to any one 
strategy in the solution of the structural design opti- 
mization problem considered. It is also shown that more 
than one local extremal design may exist in the class 
of problem treated; however, these multiple extremals 
are apparently the result of non-conve x constraint boun- 
daries. 
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In general, the multiple extremal problem may be treated 
by the warping transformation introduced by Hague in reference 
1. This approach was applied to stiffened cylinder design in 
the present study. However, the solutions converged reliably 
to a unique solution with or without the transformation. The 
exterior penalty function approach itself, when properly 
applied, is able to effectively define the global constrained 
extremal design. The multiple extremal warping transformation 
was applied to an elementary unconstrained two-variable two- 
extremal optimization problem during the study. The transfor- 
mation consistently obtained both extremal solutions. The 
optimal solutions reported here were obtained by application 
of a generalized multivariable search code, AESOP, originally 
constructed under contract to the National Aeronautics and 
Space Administration's Office of Advanced Research and Develop- 
ment. Original documentation of this code is given in 
references 1 to 3; an outline of the analysis underlying this 
code is presented below. 


MULTIVARIABLE SEARCH 

The general non-linear multivariable optimization problem 
is concerned with the maximization or minimization of a pay-off 
or performance function of the form 

<j> = ) » i = 1 » 2 N (1 ) 

Subject to an array of constraints 

c j - C-( a .) - 0, j = 1 , 2, . . . , p (2) 
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The a.j are the independent variables whose values are to be 
determined so as to maximize or minimize the performance 
function 4> ( a-j ) subject to the constraints of equation (2). 
The a.j may be looked upon as the components of a control 
vector , a, in a space of dimension N. Since maximization 
of a function is equivalent to minimization with a change 
of sign, it suffices to discuss the case in which the per- 
formance function is to be minimized. 


Multivariable optimization problems involving inequality 
constraints may also be encountered. If the constraints are 
applied directly to the independent variables 


L . . H 

a < a . < a 

i i 


(3) 


the inequality constraints define a region of the control 
space within which the solution must lie. Inequality con- 
straints on functions of the independent variables similarly 
restrict the region in which the optimal solution is to be 
obtained. In this case 


E k<«i) < MV < E i!<“i> 


(4) 


Inequality constraints can be used to restrict the search 
region directly, or, alternatively, they may be applied in an 
indirect fashion by a transformation to equality constraints. 
Several transformations may be used for this purpose. For 
example, let an equality constraint, C^, be defined by the 
transformati on 


C 


k 


< E k - e k>- E k < E k 

0 • E k < E k < E U 

' < E J - E k> a - E k < E k 


(5) 
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Constraining C ^ to zero will result in the constraint of 
equation (4) being satisfied. 


Problems involving equality constraints can be treated 
as unconstrained problems by replacing the actual performance 
function, <|>(a-j), by an augmented performance function } 
where 

P 



( 6 ) 


It can be shown that, provided the positive weighting constants 

U. are sufficiently large in magnitude, minimization of the per- 
0 

formance function subject to the cons tra i n ts , equa t i on (2-), is 
equivalent to minimization of the unconstrained penalized 
performance function defined by equation (6). This approach 
permits search techniques for finding unconstrained minima 
to be applied in the solution of constrained minima problems 
at the cost of some increased complexity in the behavior of 
the performance function, the performance response surface. 

In practical application, the weighting constants U. are 

J 

determined adaptively on the basis of response surface be- 
havior. 


Alternatives to this approach are available, notably 
Bryson's approach to the steepest-descent search, reference 
4. This method has been exploited in connection with the 
numerical solution of variational problems encountered in 
the optimization of aerospace vehicle flight paths, refer- 
ences 5, 6, and 7. However, the use of such techniques 
implies smoothness of the response surface. This smoothness 
may not be present in general; hence, the AESOP code is 
limited to the less restrictive penalty function approach of 
equation (6) 
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Numerical Solution of Non-Linear Multivariable 
Optimization Problems 

This section is devoted to a discussion of the search 
algorithms for solution of non-linear multivariable optimi- 
zation problems available in the AESOP code. A wide variety 
of search algorithms have been devised for the solution of 
multivariable optimization problems. Many of these algorithms 
are restricted to the solution of linear or quadratic problems. 
Algorithms of this type must be supplemented by more general 
search procedures if generality of solution is sought; for 
engineering problems tend to lead to non-linear formulation 
with the possibility of discontinuities in both the performance 
function response surface and its derivative. Most of the 
searches which prove effective in these problems combine a 
direction generating al gor i thm , such as steepest-descent , with 
with a one-dimensional search. Distance traversed through 
the control space in the selected direction is measured by a 
step-size, or perturbation parameter, DP. The object of the 
one-d i mens i onal search is to determine the value of DP which 
minimizes the performance function along the chosen ray and to 
establish the corresponding control vector. 

In practice, the diverse nature of non-linear multivariable 
optimization problems leads to the conclusion that no one 
search algorithm can be uniquely described as being the "best" 
in all the situations which may be encountered. Rather, a 
combination of searches, some of which may be of quite elemen- 
tary nature, provides the most reliable and economical conver- 
gence to the optimal solution. 
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One-dimensional search . Multivariable search problems 
are reduced to one-d imens i ona 1 problems whenever a search 
algorithm is used to establish a one-to-one correspondence 
between the control vector and a single scalar perturbation 
parameter, (DP). In such a situation 

a n . = a i (DP) , i = 1 , 2, . . . , N (7) 

so that equation (1) becomes 

<f> = <p(a.j) = <j> ( DP ) (8) 

Similarly, the right hand sides of equations (2) and (6) 
become functions of the scalar perturbation parameter. 

The relationship, equation (7), specifies a ray through 
the control space. As noted above, the objective of the 
one-dimensional search along this ray is to locate the value 
of DP which provides the minimum performance function value. 

Numerical search for the one-dimensional minima can be 
carried out in a local fashion, by the Newton-Raphson method, 
for example, or by a global search of the ray throughout the 
feasible region. The localized polynomial approximation is 
appropriate to the terminal convergence phase in a problem 
solution when some knowledge of the extremal's position has 
been accumulated by the preceding portion of the search and the 
problem involves a smooth function. The global search can be 
used to advantage in the opening moves of a search. In the 
early phase of a search the object is to isolate the approx- 
imate neighborhood of the minimum performance function 
value as rapidly as possible, usually with little, or no 
foreknowledge of the performance function behavior. One 
measure of the effectiveness of a search algorithm in such 
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a situation is the number of evaluations required to locate 
the minimum point to some prespecified accuracy. It can be 
shown that the most effective method of locating the minimum 
point of a general unimodal function is a Fibonacci search 3 
reference 8. In this method, the accuracy to which the mini- 
mum is to be located along the perturbation parameter axis 
must be selected prior to the commencement of the search. 

Since the accuracy required is highly dependent on the behavior 
of the performance function, this quantity is difficult to 
prespeci fy . 

Prespecification of the accuracy to which the extremal's 
position is to be located can be avoided for little loss in 
search efficiency by use of an alternative search based on the 
so-called golden section, reference 8. This is the method 
employed in the AESOP code one-dimensional search procedure. 
Search by the golden section commences with the evaluation 
of the performance function at each end of the search interval 
and at G = 2/(1 + /5) of the interval from both of these boun- 
ding points. This is illustrated in Figure 1. 

The boundary point furthest from the lowest resulting 
performance function value is discarded. The three re- 
maining points are retained, and the search continues in a 
region which is diminished in size by G. The internal point 
at which the performance function is known in the reduced 
interval will be at a distance G of the reduced interval 
from the remaining bounding point of the original interval 
for (1 - G) = G 2 . The search can, therefore, be continued 
in the reduced interval with a single additional evaluation 
of the performance function. It follows after Q evaluations 
of the performance function that the position of the extremal 
point will be known to within R of the original search region 
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FIGURE ] , SEARCH EttSED OH THE COIDKN SECTION 


where 


R = G^" 3 ) (9) 

To reduce the interval of uncertainty to .00001 of 
the original search interval, about 27 evaluations of the 
performance function are required. For a reasonable number 
of evaluations of the performance function this type of 
search is almost as efficient as a Fibonacci search. 

It should be noted that search by the golden section 
proceeds under the assumption of unimodality; hence, it will 
often fail to detect the presence of more than one minimum 
when the performance function is multimodal. If more than 
one minimum does exist, the one located depends on perfor- 
mance function behavior within the original search interval. 

Multiple Extremals on One-Dimensional Ray . The one- 
dimensional section search described above is unable to 
distinguish one local extremal from another; it will merely 
find one local extremal. This difficulty can be largely 
eliminated by the addition of some logic to the search, at 
least for moderately well behaved performance functions; 
that is, for functions having a limited number of extremals 
in the control space region of interest. An effective method 
for detecting multiple extremals is to combine the one- 
dimensional search with a random one-dimensional search on 
the same ray through the control space. This is illustrated 
in figures 2 and 3. In figure 2 the response contours of 
a performance function having two minima are illustrated 
together with the initial points used in a global one- 
dimensional search by the golden section method. The 
behavior of the function at these points is shown in 
figure 3. The left hand minimum is not apparent from 
these points. If a single random point is added in the 
interval L 0 , the probability of this point revealing the 

9 



o 



FIGURE 2 . RESPONSE SURFACE WITH TWO TROUGHS 



presence of the second minimum is 


P, * L,/L 0 (10) 

for any point in the interval AB indicates the presence of 
a local minimum somewhere in the interval AB, and any point 
in the interval BC indicates the presence of a local maximum 
somewhere in the interval BC. In this latter case, there 
must be a minimum of the function both to the left and to the 
right of the newly introduced point. 

If R random uniformly distributed points are added in the 
interval L 0 , the probability of locating the presence of 
the second minimum becomes 

P R = 1.0 - (1.0 - l 1 /l 0 ) r (11) 

The function (L-|/Lq) is a measure of the performance 
function behavior. For a given value of this behavior function 
the number of random points which must be added to the one- 
dimensional search to provide a given probability of locating 
a second minimum can be determined. 

The presence of multiple minima on a one-dimensional cut 
through an N-dimensional space does not necessarily indicate 
that the performance function possesses more than one minimum 
in a .mul ti -di mens i onal sense. It may be that the performance 
function is merely non-convex. This is illustrated by figure 
4. The performance function behavior on the one-dimensional 
search in figures 2 and 4 is identical. In figure 2 this 
indicates the presence of two local extremals; in figure 4, 
a non-convex performance function. 

When a one-dimensional search detects the presence of 
multiple extremals in the local sense above, a decision must 
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be made as to which of the apparent extremals is to be 
pursued during the remainder of the search. Here, without 
foreknowledge of the performance function behavior, logic 
must suffice. Typically, the left or right hand extremal, 
the extremal which results in the best performance, or 
even a random choice may be made. 


It should be noted that logic of this type is not cur- 
rently available in the AESOP code. The AESOP one-dimensional 
search procedure has three distinctive phases. First, each 
search algorithm defines an initial perturbation using either 
past perturbation stepsize information or a perturbation mag- 
nitude prediction as in the quadratic search below. Second, 
a perturbation stepsize doubling procedure is employed until 
a point exhibiting diminishing performance is generated. 

Third, having coarsely defined the one-dimensional extremal 
position from steps one and/or two, a golden section search 
is employed to locate the extremal with reasonable precision. 


Multiple extremals - general procedure . The multiple 
extremal search technique included in AESOP is based on 
topologically invariant warping of the performance response 
surface. The response surface is warped in a manner which 
retains all the surface extremals but alters their relative 
locations and regions of influence. The region of influence 
of an extremal is defined as the hull or collection of all 
points which lead to the extremal if a gradient path is 
followed. Reducing the region of influence of an extremal 
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FIGURE 5. FUNCTION WITH TWO EXTREMALS 



decreases the -probability of locating a point in the neigh- 
borhood of the extremal if points are chosen at random. 
Again, in an organized multivariable search, the probability 
of locating an extremal having a small region of influence 
is less than that of locating an extremal having a large 
region of influence. For example, suppose the extremals 
of the one-dimensional function of figure 5 are to be deter- 
mined in the range a L < a < by the sectioning approach. 
The four initial values employed in this technique are 
denoted by f.| to f^. 


Following evaluation at these four points, f 4 is dis- 
carded, and the function is evaluated at f 5 » Atthis point 
the right-hand extremal, e2, has been eliminated from the 
search which now inevitably proceeds to the left hand extre- 
mal at e] . 


To find the second extremal, the Function F is warped 
by writing 


F(S) = F (a) 


02 ) 
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K = (an -01 *) 


a-a* 

a^-a* 


~i 2N 


+ a* ; a > a i 


5 = -(a*-a|_) 


a*-a 

a*-a|_ 


2N 


+ a* ; a* > a (13) 


where N is a positive integer, and a* is the location of the 
1 eft hand extremal . 

A typical relationship between £ and a is shown in 
figure 6 for the case N = 1. Differentiation of equation 
(13) with respect to a when N = 1 results in 


5 ' 

V 


2 [a-a* ] 
[a H -a* ] 

2[a*-a] 
[a*-aL ] 


a a* 


a < a 


( 14 ) 


Note that as a -+ a*, £ 1 -* 0 from both the left and 
right. At a = ct|_ and at a = a H , V = 2. In the regions 
a^ < a < a* and a* < a < a^j , E, varies parabol i cal ly 
with a. Figure 6 illustrates these points. It can be 
seen that a region Aa-j centered about a* transforms into 
a smaller region A£] located in the neighborhood of 
5 = a*. On the other hand, a region Aa£ situated in the 




neighborhood of the upper search limit, , maps into a 
wider region in the neighborhood of £ = an- In general, 
the slopes at a = cq_ and a = are given by 2N ; the 
greater N, the greater the warping becomes. 

The effect of introducing a moderate warping trans- 
formation on the function of figure 5 is shown in figure 
7. It can be seen from figure 7 that the region of influ 
ence of e ] is reduced, and the region of influence of e£ 
is increased. On the warped surface search by sectioning 
commences with evaluations of performance at f' to f^. 
Following these initial evaluations f-j is discarded (as 
opposed to the discard of on the unwarped surface), 
and the function is evaluated at the additional point fg. 
The points fj and fg straddle the extremal e£ which is 
now inevitably located by further sectioning evaluations. 
Figures 8a to 8j illustrate the warping transformations 
for a range of N between 1 and 10 when the transformation 
is applied at the point a* = .5, the symmetric case. It 
can be seen that when N = 1 , twenty per cent of the 
warped control space corresponds to approximately 45 per 
cent of the unwarped control space in the vicinity of 
the transformation origin (cx = .5). When N = 10 twenty 
per cent of the warped control space transforms into 
ninety per cent of the unwarped control space. 


Sectioning Parallel to the Axes . The independent 
variable perturbation algorithm in the sectioning search 
i s 

Aa.j = 0, i / r 

=DP, i=r r = 1 , 2 , . . . , N (15) 
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This is simply the parametric or univariate search 
approach. All but one of the independent variables are held 
constant while a one-dimensional search parallel to the Rth 
variable axis determines the best value of the remaining 
variable, a r . The variable a r is then set to this value t 
and the process is repeated with one of the remaining 
independent variables. When all N independent variables 
have been perturbed in this way, a sectioning search cycle 
has been completed. 

The N-dimensi onal search can then be continued with 
another cycle of sectioning or by one of the other search 
techniques described below. In practice, it has been found 
advantageous to -perturb the independent variables in a 
random order within each sectioning cycle. The method can 
be used in conjunction with either a local or a global 
search as outlined in the two preceding sections. The 
behavior of this search in the solution of a straightforward 
two-variable optimization problem is illustrated in Figure 9. 
It may be noted that the AESOP code searches from boundary 
to boundary in each variable using a golden section search 
procedure . 


Sectioning to Define Local Sensitivities . The sec- 
tioning search can readily be applied to the problem of 
performance or constraint sensitivity determination. Thus, 
by the device of omitting the updating of each control 
variable a r following the sectioning search on the r 1 -* 1 
parameter, the sequence of sectioning searches is performed 
about a fixed nominal point. When such a search is per- 
formed in the vicinity of a known extremal point, the 
penalties for off-optimal design can be assessed. Away 
from an extremal point, the search merely provides local 
sensitivities in a similar manner to the manual perturbation 
methods employed in conventional trial and error design . ' “T~~ 
evolution. 23 jC?* ’ 

■j 
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The steepest-descent search 


Steepest-Descent Search , 
algorithm is 

u<*> -:wr'j{|i} - -‘{k 2 }J 

x /(DP) 2 - DC CK3] -1 {DCT 
Kt K 2 [K 3 ]' 1 {K2> 

- [W]- 1 [§|] T [K 3 ] j {DC} (16) 


Here, the matrix W is the metria tensor of the control space 
and serves to define a generalized measure for the magnitude 
of a control vector perturbation. The vectors { 9<f>/ 3ct> and 
{3C/3a} are defined as 

3 <p 3 ()> 3 <p 

3a! 5 3a 2 ’ ‘ * * 3« n 


and 


3 C 


j 9ot l ’ 

respecti vely . 


3 C 3 C 

3a, ’ ‘ ‘ 3a„i 

c n I 

The K matrices are defined as 


K i ■ Lfil ‘If } < ,7) 

{ K 2 }= (18) 

[K 3 ]= [fllCW]- 1 [||l T (19) 
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The perturbation parameter, (DP), is defined by 


(DP) = LH [WJtAcx} (20) 

The vector is the desired change in the constraint 
functions. For an unconstrained problem, (16) reduces to 

{Aa} = -[W]- 1 ^} Aj^p- 2 (21) 


The performance function change associated with the pertur 
bat ion of equation (16) is 



LK2j[K 3 3 


')■'(' 


{K ? }V (DP) 2 - IDC) [K 3 ] 


-1 



+ L K 2j t k 3 3 1 (22) 


Equation (16) does not specify a one-dimensional 
search directly since the perturbation parameter (DP) 
and each component of the constraint vector change DC 
can be independently specified. This difficulty is 
conveniently eliminated if the components of DC are 
expressed in terms of the perturbation parameter. Let 
(DP) and DC be arbitrarily assigned, say (DPq) and DCq , 
respectively. Now consider the one parameter set of 
values for DC defined by 

= (dPq^ * DC 0 ( 23 ) 

It follows from equations (16) and (22) that (23) spec- 
ifies a one parameter family of perturbations in which 
the non-linear performance and constraint functions vary 
linearly with (DP) t to the first order. 
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Equations (16) to (22) are valid for small pertur- 
bations in the independent variables provided the der- 
ivatives involved are continuous in the region of the 
control space defined by equation (20). In practice, 
when this condition is not satisfied, the steepest- 
descent algorithm can be used to locate a promising 
direction for a one-dimensional search provided the 
derivatives are computed numerically. In this case, 
however, equation (22) ceases to provide an accurate 
indication of performance function behavior along the 
specified ray. 

When dealing with performance and constraint 
functions having continuous first derivatives, the 
perturbation parameter value to be used in equation (16) 
can be determined from a second order Taylor's expan- 
sion of the performance function behavior in terms of 
DP. The coefficients in this series expansion can be 
readily obtained from the conditions of zero change for 
DP = 0, linear slope for DP = 0, and from the actual value 
of the performance function at a point in the neighborhood 
of the point at DP = 0. This method for determining the 
best perturbation parameter value is discussed in some 
detail in references 5 and 6. When dealing with less reg- 
ular functions, the one-dimensional search by sectioning 
can be used to determine the perturbation parameter value. 
This is the technique employed in the optimization program, 
AESOP, references 1 and 2; for the AESOP code converts all 



constrained optimization problems to unconstrained problems 
by the penalty function device, equation (6). The resulting 
response surface combines both performance function and 
weighted constraint functions. Inevitably, this surface has 
a more complex topology than that of the unconstrained per- 
formance function. Program AESOP is also limited to the 
penalty function approach to constrained optimization, and, 
hence, it utilizes the reduced algorithm of equation (21) 
rather than the explicit constraint algorithm of equation 16. 


Steepest-Descent Weighting Matrices . The weighting 
matrix introduced in equations (16) and (20) must be positive 
definite to assure a positive distance between any two non- 
coincident points in the control space. Apart from this 
restriction, the choice of weighting matrix is arbitrary. 
Inspection of equation (16) reveals that any descending 
direction is a steepest-descent path for some choice of 
the weighting matrix W. This can be simply illustrated 
when only two independent variables are involved. Figure 10 
depicts a small region of the control space R 2 . The per- 
formance function response contours appear as a series of 
parallel lines on this microscopic region of the control 
space. The perturbation zones corresponding to three 
weighting matrix choices are shown. The first zone corres- 
ponds to the choice of a unit matrix for W. It follows from 
equation (20) that for a given value of (DP) 2 the search 
zone is a circle of radius (DP). The steepest-descent 
direction is that in which the performance improvement is 
greatest. This is the direction of a line from the origin 
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of the circular search zone to that point on its circumference 
which provides the smallest value of the performance function 
<t>(a). With this choice of weighting matrix, the steepest- 
descent direction is perpendicular to the response contours. 
Paths of this type are illustrated in figure 11 by the solid 
lines emanating from points A and B. From the nominal point 
A, search perpendicular to the performance response contours 
is very efficient. From r-oint B, however, this type of 
search results in the meandering path illustrated. It is 
assumed here that once a steepest-descent direction is located, 
an exhaustive search for the minimum in that direction will 
be undertaken in view of the high cost of recomputing the 
derivatives in many problems. Even if this were not the 
case, search normal to the response contours can often be 
improved upon. For example, it is obvious that even in the 
straightforward two-dimensional problem of figure 11 the 
dashed search direction is superior. This direction requires 
a priori knowledge of the extremal's position, information 
not normally available. 


Returning to figure 10, the second search zone 
depicted corresponds to the choice of a diagonal matrix 
for W. The positive-definite constraint on W requires 
that all diagonal elements of the weighting matrix be 
positive. In this case the search zone becomes elliptical 
with the major and minor axes of the ellipse being parallel 
to the coordinate axes. It may be noted that as either of 
the diagonal elements of W becomes large in relation to the 
remaining element, the corresponding element in W inverse 
together with the predicted change in the associated inde- 
pendent variable becomes small. In the limit this reduces 
the search to a one-dimensional search in the remaining 
coordinate. The perturbation zone then becomes a slit 
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FIGURE IQ, -PERTURBATION ZONES CORRESPONDING 
TO THREE WEIGHTING MATRICES 



parallel to that coordinate axis of length 2 • (DP), as 
illustrated in figure 10. In the case illustrated, the 
steepest-descent path is in the descending a-j direction. 

Finally, the search zone correspondi ng to the choice 
of an arbitrary positive-definite weighting matrix is 
shown. From equation (20) and the positive-definite con- 
straint on W, the search zone remains elliptical, but the 
principle axes may now have an arbitrary orientation to 
the axes of a-j and ag. It follows that since the elliptic 
search zone can have any orientation and eccentricity, any 
direction in the control space is a possible steepest- 
descent path j for in all cases, the path of steepest-descent 
lies in the direction of a line joining the search zone 
origin to the lower point of tangency between the boundary 
of the search zone and the performance function response 
contours. The discussion above may readily be extended 
to control spaces of higher dimensionality. 

When attempting the solution of optimization problems 
by the steepest-descent method, the analyst is constantly 
faced with the problem of choosing a satisfactory weighting 
matrix for the search continuation. The problem is com- 
pounded by the fact that the slopes of the performance 
function with respect to the independent variables can, 
and frequently do, vary by many orders of magnitude. The 
arbitrary choice of a unit matrix in such situations can 
lead to distressingly slow convergence of the numerical 
search; for it is in the nature of many problems that in 
those directions in which the slopes are greatest the re- 
sponse surface is highly non-linear. Only small pertur- 
bations will be successful in the direction of these strong 
control variables. In those directions in which the slopes 
are small, the contours are often relatively linear, and 
large perturbations may be required in these weak control 

variables . 
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In such situations the local steepest-descent direction 
for [W] = [I] is quite misleading; for contrary to the resul- 
ting steepest-descent direction which, by equation (21) 
results in independent variable perturbations which are in 
proportion to the response surface partial derivatives, the 
best direction in which to proceed may well involve large 
perturbations in the weak control variables of small slope. 
This behavior is illustrated for a two-dimensional case in 
figure 11 by the dashed line emanating from B. 

The problem of choosing a satisfactory weighting matrix 
also arises when the steepest-descent search is applied in 
its variational form, reference 5, and when a combination of 
continuous control variables and parameters are encountered 
as in the optimization of multiple-arc problems in flight 
path optimization problems, reference 6. In these references 
it is suggested that the weighting matrices be based on the 
first derivatives of the unconstrained performance function 
with respect to the control. This approach can be used in 
the solution of multivariable optimization problems also, by 
writing 



= 0, i f j 


In practice, alternate use of the resulting combined weighting 
matrix and the unit matrix tends to provide a reasonable con- 
vergence rate at points well removed from the extremal. The 
AESOP code employs such a matrix in combination with a search 
range non-dimensi onal i zati on term and a learning factor. This 
factor emphasizes perturbation of control parameters which 
change in a mo no tonic direction and de- emphasizes those per- 
turbations whose perturbations fluctuate in sign. 

Random Ray Search . The difficulty, in some cases, of de- 
fining a suitable control variable metric tensor together with 
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the fact that any descending path is a steepest-descent 
direction for some choice of metric tensor suggests the 
possibility of searching along a random ray through the 
control space. The algorithm for random ray search is 

ha. = R. (+DP) , i = 1 , 2, . . . ,N (25) 

where the , proporti onal to the direction cosines of the 
ray, are uniformly distributed random numbers satisfying 


-1 .0 < R. < +1 .0, i = 1 , 2, . . . , N 


The positive sign in equation (25) is taken if ^ pp y is 
negative; the negative sign is taken when this derivative 
is positive. 


The utility of this type of search tends to be in 
proportion to the complexity of the performance function 
response contours. On a well-behaved problem there is 
little to recommend this type of search; on a problem 
involving unexpected behavior on the part of the performance 
function, a random ray search can be quite efficient, par- 
ticularly when used with the pattern search acceleration 
procedure below. The method is, of course, equivalent to 
a steepest-descent search using a randomly generated metric 


tensor . 


Quadratic search . An alternative systematic approach 
to the definition of an arbitrary or empirical weighting 
matrix is provided by second order or quadratic method. It 
can be shown, for example, in reference 1, that on an 
elliptic second order response surface the weighting 
matri x 


W ij 




( 26 ) 


will immediately define the point 


* Also known as the Newton-Raphson method 
34 



{a*} = {a Q } + {& a} 


(27) 


where {6a} is computed from equation (21) with (DP) 2 = .5K], 
On a more general non-linear response surface, equation (27) 
merely defines a direction for subsequent search in the 
manner of the steepest-descent technique. This is illus- 
trated in figure 12. Here, the approximating elliptical 
contours computed at point 0 define an approximate extremal 
location at P through equations (26) and (27). Subsequent 
search along the ray OP results in the definition of a one- 
dimensional extremal. This point is then used to fit another 
approximating elliptic contour, and the process is repeated 
until the extremal point at Q is located. 

The quadratic search procedure can be quite rapid in 
control spaces of low dimensionality. In high order spaces 
the approach is usually impractical as a result of the 
requirement to establish the second order weighting matrix 
of equation (26). In many practical engineering problems 
these derivatives cannot be obtained in closed form; in such 
cases the derivatives must be obtained numerically, for 
example, reference 1. Computation of these derivatives 
requires at least (N + l)(N+2)/2 evaluations of <J> at each 
point where an approximating quadratic is employed. Clearly, 
for large N this computation may become impractical. 

Davidon or FI etcher-Powel 1 Method . Davidon's method 
is a hybrid first order/second order technique. The objec- 
tive of Davidon's method is to arrive at a reasonable 
approximation to the second order weighting matrix of 
equation (26) without the use of (N+l)(N+2)/2 evaluations 
of . It can be shown that on a quadratic (second order) 
response surface N steepest-descent searches performed in 
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the manner described previously will lead to definition of 
the weighting matrix of equation (26), if the following 
formula is employed: 


[W]^! = [W]/ + [A] 1 + [B]. 


(28) 


where 



(29) 


(30) 


[Wlf 1 = [I] (31) 

Here, |_Acx J i is the change in position during the i^* 1 one- 

dimensional search and 


A . i&l 

L 3a_|. 


is the change in gradient vector between the beginning and 

J_ L. 

end of the i l one-dimensional search. On a numerically 
well-behaved function this technique may work well. When 
appreciable numerical noise is present in the calculation, 
the method may produce erratic convergence to the extremal 
point, or convergence failure. Again, {1^-} is the gradient of 
* defined as {f£_, ||g • • • • 

Pattern Search . In the present report, pattern search 
refers to a search which exploits a gross direction revealed 
by one of the other searches. The search algorithm is 


Aa . = ( a ? - a* ) 
i v i 


(DP), 


i - 1, 2, 


N (32) 


V; 


n 

/#? 


where a? and a| are the components of the control vector 
before and after the use of a preceding search technique. 
This type of search is illustrated in figure 13 following 
a section search. The combination of a section search 
and a pattern search in the problem illustrated leads 
directly to the neighborhood of the extremal. Repeated 
sectioning, on the other hand, would be a very slowly 
converging process due to the orientation of the contours 
with respect to the axes of the independent variables. 

It may be noted that a simple rotation of the independent 
variable axes by 45° results in sectioning alone becoming 
a rapidly converging process in this example. The pattern 
search can also be used to accelerate the steepest-descent 
process provided it follows two successive descents as in 
figure 14. 

Adaptive Search . Adaptive search is a form of small 
scale sectioning; however, instead of locating the position 
of the one-dimensional extremal on each section parallel 
to a coordinate axis, the coordinate is merely perturbed 
by a small amount, Aa r , in the descending direction. 

The search commences with a small perturbation in one 
of the independent variables, a r ; a positive perturbation 
is first made; if this fails to produce a performance 

improvement, then a negative perturbation is tried. If 
neither of the perturbations produces an improved perfor- 
mance value, the variable retains its nominal value, and 
A<* r is halved. If a favorable perturbation is found, the 
variable a r is set to this value, and Aa r is doubled. The 
process is repeated for each independent variable in turn, 
the order in which the variables are perturbed being 
chosen randomly . At this point an adaptive search cycle 
is complete, and the cycle is then repeated. A two- 



FIGURE 13. PATTERN SEARCH FOLLOWING SECTIONING 
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dimensional illustration of this search is presented in 
figure 15. In the particular problem illustrated, the 
method converges rapidly reaching the neighborhood of the 
extremal within six evaluations. 


The search algorithm can be written in the form 


Aa r = 2.0^ S r' T ^^ . (DP) 


(33) 


where S r is the number of cycles in which the search has 
successfully perturbed the r^h independent variable, and 
T r is the number of cycles in which a perturbation of the 
r't' 1 variable has proved unsuccessful. While this search 
can be looked upon as a one-dimensional approach, this 
viewpoint is somewhat artificial. Here, the scalar quantity 
(DP) merely defines an initial perturbation for each inde- 
pendent variable. Once started the search proceeds inevi- 
tably to its conclusion, the perturbation in each independent 
variable being adaptively determined according to equation 
(33) on the basis of the performance function response 
contour behavior encountered during the particular problem 
solution. This search can be quite efficient when used in 
combination with the pattern search acceleration procedure. 

Magni f i cati on . When studying discrete models of con- 
tinuous systems of the type encountered in certain engineering 
problems such as aerodynamic shaping or structural design 
problems, there is a tendency on the part of some search 
algorithms to achieve a favorable shape before satisfying 
the desired constraint levels. In such cases, when it is 
known that the unconstrained extremal is the null vector, 
a simple magnification search can lead to rapid convergence 
to the desired solution. The magnification algorithm is 
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Here (DP) is positive and all components of the control 
vector are to be simultaneously perturbed. Generally, 
the unconstrained extremal point corresponds to the null 
vector; this method may prove efficient. 

Arbitrary Ray Search . In practical design optimization 
a search along an arbitrary multidimensional ray can be of 
utility. For example, when two minimal extremal solutions 
appear to be possible, a search on the ray connecting the 
two points should reveal the presence of a maximal extremal 
somewhere on the ray between the two minimal extremals. The 
algorithm for this search is 

Aa. = (a i - a.j ) (DP), i = 1, 2, . . ., N (35) 

where and a\ are the two minimal extremal points. In 
* 1 *2 

general, a., and a., may be any two points in the control space. 

Random Point Search . A straightforward Monte-Carlo 
search which examines point designs distributed in a uniform 
random manner within the feasible region is often of utility 
when the response surface is of a complex nature. Such a 
search is included in the AESOP code primarily for use as a 
nominal point design generation procedure. 

STRUCTURAL ANALYSIS OF A STIFFENED CYLINDER 

The structural analysis employed in this study is iden- 
tical to that employed in a previous National Aeronautics 
and Space Administration-sponsored study reported by Morrow 
and Schmit in reference 9. Morrow and Schmit's analysis of 
a stiffened cylinder is reproduced for completeness in Appen- 
dices A and B. The computer code for the failure mode analysis 
employed in the present study is that developed in reference 
9 without change. However, the FI etcher-Powel 1 / F i acco- 
McCormick procedure of reference 9 was replaced by the refer- 
ences 1 through 3 optimization program AESOP during this study. 



Figure 16 illustrates the class of structures considered. 
The structure is a cylinder stiffened in both longitudinal 
and circumferential directions. A typical shell element is 
presented in figure 17. The shell element is defined by 
seven parameters: 

d x = Depth of longitudinal stiffeners * 
d^ = Depth of ci rcumferenti al stiffeners* 
j i = Spacing of circumferential stiffeners 

A 

SL^ = Spacing of longitudinal stiffeners 

t = Skin thickness 
s 

t = Thickness of longitudinal stiffeners 

t, = Thickness of circumferential stiffeners 
9 

These parameters are the design vector components in the 
stiffened cylinder design optimization problem considered. 

L a J = LV t x* S* d x* d 4>’ A x* ( 36 ) 

The optimization problem considered is that of weight 
minimization subject to yield and buckling failure mode 
constraints. Cylinder weight, the payoff function , is 
taken from Appendix C of reference 9. 


* Positive number means inside stiffeners; negative number 
means outside stiffeners. 
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(37) 

Multiple load cases are considered. For each load case, 
inequality constraints are placed on five yield failure 
modes. Equations A27 and A29. 

1 . Skin yield 

2. Longitudinal stiffener yield in tension 

3. Longitudinal stiffener yield in compression 

4. Circumferential stiffener yield in tension 

5. Circumferential stiffener yield in compression 

A distortion energy type criterion is used in the skin for 
the biaxial state of stress. In the stiffeners, the uni- 
axial state of stress must have a value between the compres- 
sion yield value and the tension yield value. 

Three buckling failure modes of the cylindrical shell 
are considered: 

1. Gross buckling of the entire cylinder 
(gross buckling) 

2. Buckling of the cylinder between circum- 
ferential stiffeners (panel buckling) 

3. Buckling of the cylindrical skin (skin 
buckl i ng) 

In the gross buckling mode, the effect of the stiffeners 
are averaged over the stiffener spacing. For the panel 
buckling mode, only the longitudinal stiffeners are 
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averaged in. Bending stiffnesses of the sti f f eners , the 
torsional stiffnesses, and the effects of eccentricity are 
taken into account. The cylinder buckling analysis is a 
linear classical small displacement analysis, assuming 
simply supported boundaries and a uniform prebuckled mem- 
brane force and displacement distribution. The same anal- 
ysis is used to determine the critical loads for gross, 
panel, and skin buckling, by substituting the appropriate 
stiffness properties and displacement patterns. An expres- 
sion for the buckling load in terms of the mode shape is 
given by equation A20 or A22. The critical buckling load 
is found by determining the buckling loads for a large 
number of mode shapes and selecting the lowest of these 
loads as the critical value. 

Stresses and strains in the skin and stiffeners prior 
to buckling are determined from the membrane force distri- 
busion, equations A 1 2 , A13, and A14. The strains in the 
stiffeners where they join the skin are assumed to be the 
same as the corresponding strains in the skin. 

Three failure modes are considered for the stiffeners: 

1. The longitudinal buckling stress is 
calculated from equation A23 

2. Outside circumferential stiffeners can 
buckle either when the cylinder expands 
or contracts under load 

3. Inside circumferential stiffeners can 
buckle only when the cylinder contracts. 
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The expression for the circumferential stiffener 
critical strain, equation A24, derived in Appendix A, 

Section A. 9, is verified for two limiting cases in 
Appendix B. In the stiffener buckling analysis, simply 
supported boundaries are assumed at all edges where the 
stiffener connects with the shell or the other stiffeners. 

In the optimization procedure employed herein, all 
buckling and yield constraints are expressed in the form 

P k < P k (38) 

where P|< is the load which the cylinder is to support, and 
is the critical load for the k*' 1 failure criteria. These 
inequality constraints are converted to equality constraints 
by a one-sided transformation in an analogous manner to 
that of equation (5). 


0; P k < P k 

(39) 

(P k - P k ) 2 i P k > P k 

(40) 


Equations (30) and (31) may define a large number of con- 
straints as k varies over all failure criteria and load 
conditions. For example, with fifty longitudinal and fifty 
circumferential buckling modes considered for gross, panel 
and skin buckling and ten load conditions, the number of 
constraints to be considered is 

N k - N fc x Njj x n l (41) 


where 


= Number of constraints 
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Np L = Number of failure criteria considered 
equals three 

N_.j = Number of buckling mode combinations 
equals 50 x 50 = 2,500 

= Number of load conditions equals ten 

In this particular example, there would be 75000 buckling 
constraints. Each constraint corresponds to a particular 
mode shape and load condition. In applications to-date, 
the maximum number of load conditions considered is three 
which, in the example, would result in the maximum of 22500 
constraints . 

In practice, computation of the violations for this 
number of constraints becomes a time consuming process 
when each constraint violation is checked at a number of 
points in the control space. In reference 9 Morrow and 
Schmit reduced the computational time by introducing partial 
or approximate analyses. At the beginning of a synthesis 
a complete cylinder analysis is performed. The word 
"complete" means that a large number of buckling mode 
shapes are examined. (Each mode shape is a constraint). 

The subset of mode shapes that is most active in the com- 
plete analysis is saved, and these mode shapes are examined 
in several subsequent approximate analyses. These approx- 
imate analyses are carried out during a succession of 
moves through design variable space. Periodically, a 
complete analysis is performed, and the subset of mode 
shapes used in the approximate analyses is redefined. An 
approximate analysis is approximate only in that the number 
of buckling mode shapes examined -j s small compared with the 
number of buckling mode shapes examined on a complete analysis. 


50 


This same procedure was employed in the present study. 
Following Morrow and Schmit, when a large number of constraints 
are checked, the analysis is defined to be complete; otherwise 
the analysis is approximate. In a typical load case problem a 
complete cylinder analysis required three seconds on a CDC 
6600 computer while an approximate analysis in which 70 
active constraints were retained required .03 seconds. 

Clearly, computer time requi rements can be significantly 
reduced by a good mix of complete and approximate analyses. 
However, it must be noted that when successive complete 
analyses are too widely separated in the control space, 
convergence rates may diminish to the point where compu- 
tational time gains are negated. 


MINIMUM WEIGHT DESIGN OF 
STIFFENED CYLINDERS 

In this section a variety of minimum weight stiffened 
cylinder designs are considered. Cylinder designs are 
optimized by a coupled version of the multivariable 
optimization program AESOP of references 1 through 3 and 
the Morrow and Schmit structural analysis program of 
reference 9. The technique employed is that of design 
evolution by repetitive perturbation and analysis as 
illustrated in figure 18. Here, a nominal design charac- 
terized by the control vector cr„ is supplied to the 
optimization program. The design parameters corresponding 
to this control vector are, in turn, passed to the system 
model which is, in this case, a stiffened cylinder analy- 
sis program, by the optimizer. The system model operates 
in a black box fashion and evaluates the system performance 
characteristics consisting of the weight and failure 
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FIGURE 18. OPTIMIZER SCHEMATIC 
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criteria constraints. No explicit derivative computations 
are required by the optimizer. The performance character- 
istics are then transferred from the system model back to 
the optimization program. 

At this point a succession of perturbed designs are 
defined within the optimization program, and the computa- 
tion loop is repeated with improving designs being 
retained and inferior designs being rejected — an evolu- 
tionary process. The optimization program contains a 
variety of perturbation algorithms as illustrated in 
figure 19. These algorithms may be employed separately 
or in combination at the analyst's option. The repetitive 
design process outlined defines a sequence of gradually 
improving designs which lead from an arbitrarily selected 
initial design which may or may not satisfy the constraints 
to a locally minimal weight design which satisfies the 
constraints. The initial design need not satisfy the con- 
straints because an exterior penalty function is used in AESOP. 

Search Sequence 

In general, a combination of search algorithms will 
tend to produce more reliable convergence in the solution 
of a non-linear parameter optimization problem than the 
repetitive applications of any single algorithm. Through- 
out the remainder of this report the search algorithm 
combination employed will be designated by an array, M. 

The algorithm corresponding to an element of M, say M-j , 
is obtained from the following set of values: 

M.j = 1 , Sectioning Search 
= 2, Pattern Search 
= 3, Magnification Search 
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FIGURE 19. 


ALGORITHMS FOR OPTIMIZATION 














= 4, Steepest-Descent Search 
= 5, Adaptive Creeping Search 
= 6, Quadratic (Newton-Raphson) Search 
= 7, Davidon (FI etcher-Powel 1 ) Search 
= 8, Random Point (Monte Carlo) Search 
= 9, Random Ray Search 


Typical Optimization Algorithm Behavior 
in Stiffened Cylinder Design 

The relative efficiencies of a variety of optimization 
algorithms was studied by repeated solution of one stiffened 
cylinder design problem. An aluminum design was used in 
the study. Cylinder physical characteristics were: 

L = 291, Cylinder length, inches 
R = 95.5, Cylinder radius, inches 
E = 10.5 x 1 0 6 , Young's modulus of elasticity 
1 b s / i n 2 

Y = .101, Weight density, lbs/ in 3 

v = .33, Poisson's ratio 

Oy = 50,000., Yield stress, lbs/in 2 

A single load condition was considered. This was: 

N = 800., Applied axial force per unit length 
of circumference, lbs/ in. 

P = 0.0, Radial pressure 

The problem is solved both with and without constraints 
on the stiffener depth/ thi ckness (d/t) ratios. When the 
(d/t) limits are omitted t the problem considered is iden- 
tical to that of problem 7-1 in reference 9 . 
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Solutions Without (d/t) Limits 


Table I (figure 20) presents the results obtained 
when (d/t) limits are omitted. In the lower portion of 
the table, solutions obtained by Morrow and Schmit in 
reference 9 are included for comparison purposes. The 
table presents both final design vector and the corres- 
ponding weight obtained. The solutions obtained by 
program AESOP all utilize the starting vector used in 
Case 7-1 of reference 9, that is: 

La qJ = ^ » .05, 1., 2., 8., 3 .J (42) 

Starting vector elements are defined in equation (27). 
Cylinder weight corresponding to the starting vector is 
1681.7 lbs. 

It can be seen from figure 20 that a variety of mini- 
mal cylinder weight designs are obtained from both the 
AESOP search algorithms and the method of reference 9. 

Two search combinations achieve weights well under the 
solutions reported in reference 9. These combinations 
are 

|M.J = |_9, 2, 9, 2, 5, 2, 3j 

and 

LM.J = L 7 > 2, 7, 2j 


The first combination is the recommended procedure for 
the solution of stiffened cylinder problems. These two 
search algorithms produce minimum weight designs of 682 
and 684 lbs., respectively, approximately SO % lighter than 
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the best corresponding solution in reference 9. Depth/ 
thickness ratios utilized in the two best AESOP designs 
are 

(d/t) ^ = 10 6 and 320 
(d/t) x = 14.1 and 13.9 

It can be seen that in both cases completely unrealistic 
circumferential stiffener (d/t) values are being employed. 

The designs attained appear to be quite different from the 
heavier cylinder of reference 9. Detailed results of solutions 
without d/t limits are presented in Tables D1-D7 of Appendix D. 

Solution With (d/t) Limits Imposed 

Unrealistic stiffener depth-to-thi ckness ratios can 
be eliminated by the introduction of inequality constraints 
on stiffener section geometries. In the present report 
inequalities of this type are transformed into equality 
constraints by constraining F x and F^ to zero, where 

F x = 0; (d/t) x « 20 

= C(d/t) x - 20] 2 ; (d/t) x > 20 (43) 

- 0; (d/t) t « 20 

= [(d/t) + - 20] 2 ; (d/t)^ > 20 (44) 
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Figure 21, Table II, presents solutions obtained with 
these constraints imposed. The constraints are imposed 
as "soft" boundaries by the device of limiting the 
growth of the adaptively determined constraint weighting 
factors, (U. of equation 6), within the optimizing pro- 

J 

gram AESOP. Detailed results of solutions with d/t limits 
imposed are presented in Tables D8-D13 of Appendix D. 

Introduction of (d/t) limits reduces the spread of 
minimal designs, figure 21, which now exhibit a range 
of weights ranging from 807 to 894 pounds. It should 
be noted that all designs in figure 21 produce lower 
weight than the reference 9 designs of figure 20 which 
are not subject to the additional (d/t) constraints. 

The weight improvement over these designs varies from 
eight per cent to seventeen per cent when compared to 
the best reference 9 design , (979 pounds). 

For all AESOP solutions shown in figures 20 and 21, 
except for the quadratic search in figure 20, two suc- 
cessive runs of approximately ninety system seconds each 
were made. Therefore, the results obtained indicate how 
well the various combinations of search procedures did 
for a given amount of computer time. 

The spread of weights shown do not indicate that 
AESOP has obtained multiple solutions but instead indicate 
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the design vector obtained after a given amount of com- 
puter time. It is obvious that the AESOP solutions 
shown in figure 20 have not all converged. The quadratic 
search without d/t limits terminated during the first 
run due to a singular matrix. 

It can be seen from figure 21 that after imposing 
constraints on the stiffener depth to thickness ratio, 
the selection of search procedure used is less critical. 


Nature of the Multiple Extremal 
Minimum Weight Designs 

Physically, the existence of multiple minimal 
weight designs in the design space is unlikely unless 
the multiple extremals are constraint-induced; for 
the cylinder weight diminishes monotonically with stif- 
fener thickness, stiffener depth, and stiffener spacing. 
Now the constraint boundaries to the stiffened cylinder 
problem are probably highly non-convex since they in- 
clude multiple buckling criteria. (One can observe a 
simple boundary of this type for conical shells on page 
511 of reference 15 ). 
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Figure 22 illustrates one possible form of a non- 
convex constraint boundary problem in the case of a 
hypothetical two-dimensional problem. Here, the perfor- 
mance contours vary smoothly, but the constraint boundaries 
involve a sequence of non-convex arcs. The forbidden 
region is determined by the shaded region. In the problem 
illustrated, the constraints introduce five possible mini- 
mal solutions. The global extremal is the extreme left 
minimal solution. Two typical exterior penalty function 
convergence paths are shown in figure 22. In both paths, 
the constraint boundaries are pierced early in the search. 
At this stage in the optimization process, the exterior 
penalty function procedure is concentrating on payoff 
function improvement , and the constraint violation is of 
secondary importance. In the lower path, pursuit of 
improving performance leads to Point A before increasing 
weight on constraint violations results in a gradual 
loss of performance to satisfy the active constraints. 

This path leads to the global optimal at Point B. A 
significant point regarding constraint violation is that 
the constraint initially violated at C,and two other inter- 
mediately active constraints along arc CA, are no longer 
active at Point A. A second typical exterior penalty 
function path is illustrated by arc CD. Here, the con- 
straint weights have been increased more rapidly than on 
arc CA, and an inferior extremal is located at Point D. 
Relaxation of the constraint weights at D will result in 
further payoff function improvement along a path such as 
DB. 
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On a complex response surface or in a search from a 
badly chosen nominal, convergence to the global extremal 
may require several constraint weight relaxations. As a 
matter of practice during this study the constraint weights 
were reduced at least once on each problem. 

It should be noted that when high constraint weights 
are employed, the exterior penalty function method will 
usually fail to find the global extremal of a constrained 
optimization problem such as that illustrated in figure 22. 
With high constraint weights, the search turns to follow 
a constraint boundary as it is encountered. In conse- 
quence, the nearest constrained extremal point is found; 
this may or may not be the global extremal. Similar results 
may be encountered when an interior penalty function tech- 
nique is employed with low constraint weights. This may 
well be the reason for some of the very high weight optimal 
cylinder designs obtained in reference 9. 

When the unconstrained payoff function response sur- 
face itself possesses more than one extremal, successive 
constraint weight reductions may or may not lead to the 
global extremal. In such cases, the multiple extremal 
procedure of reference 2 may be combined with the constraint 
relaxation procedure in a search for successive minima. 

It is apparent that the multiple "extremal" solutions 
of figure 20 could be attributed to non-convex buckling 
constraint boundaries. It will be shown in a later section 
that these constraint boundaries do have the general 
characteristics exhibited in figure 22 and that the vari- 
ation in unconstrained performance along a line such as 
DB in figure 22 is essentially a smooth monotonic improve- 
ment to the minimal weight design point. 
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FIGURE 22. CONSTRAINT INDUCED EXTREMALS 


Comparison Between Interior and Exterior 
Penalty Function Approach to 
Stiffened Cylinder Design 

In this section a series of direct comparisons between 
minimum weight cylinder designs obtained by an interior and 
an exterior penalty function optimization approach are 
presented. The interior penalty function method is that of 
Fiacco and McCormick applied through the FI etcher-Powel 1 
unconstrained minimization algorithm as reported in refer- 
ence 9. The exterior penalty function method employed is 
that of program AESOP, references 1 through 3, which may be 
utilized with any of the search algorithms of figure 19. 

When using the exterior penalty function approach, the con- 
straint weights were routinely relaxed onae in each solution 
as discussed in the previous section. In one particularly 
difficult problem, the constraint weights were relaxed four 
times. This problem is discussed in some detail below. All 
exterior penalty function solutions obtained in this section 
used a combination of searches. 

[M.J =[9, 2, 9, 2, 5, 2, 3j 

Case identification conforms to that of reference 9. 

Each of the interior penalty function cylinder designs in 
reference 9 is considered; however, not all of the starting 
solutions of reference g are employed in the exterior pen- 
alty function solutions. Details of the results presented 
in this section are given in Appendix C. 
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Case 1-1, Three Load Cases, Cylinder Length = 165" , 
Cylinder Radius = 60" . Load conditions for this problem 
are : 


1. 

Axial load 

700 

1 bs/in. , 

External 

pressure 

o 

« 

o 

2. 

Axial load 
lbs/in 2 

940 

1 bs/in. , 

External 

pressure 

-2. 

3. 

Axial load 
1 bs/i n 2 

212 

1 bs/in. , 

External 

pressure 

0.4 


A summary of the results obtained for Case 1-1 is presented 
in Appendix C, Table Cl. As noted, two searches were made 
on the computer to obtain convergence with a constraint 
weight reduction between the two searches. The first search 
consisted of 15 optimization cycles. (A cycle consists of 
one application of the search algorithms specified). The 
second search consisted of ten optimization cycles. The 
search techniques used in each optimization cycle were 
Random Ray, Pattern, Adaptive Creeping, and Magnification 
in the combination defined above. Values of 10.0 were 
arbitrarily selected as initial constraint weights for both 
searches . 

The final values of the design variables obtained by 
the AESOP exterior penalty function, figure 23, and the 
corresponding cylinder weight 226 lbs. indicates that the 
AESOP solution has converged to essentially the same con- 
figuration as shown on page 71 of reference 9 , (cylinder 
weight of 231 lbs.). Table C2, Appendix C, presents the 
detailed results obtained by AESOP for Case 1-1 when the 
final values of the design variables shown on page 71 of 
reference 9 are used as the nominal values for AESOP. 

It can be seen that a small performance gain is possible. 
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Table C3 presents the results obtained by AESOP for 
Case 1-1 when the absolute values of all behavior vari- 
ables were constrained to be less than or equal to 1.0 
This procedure was followed in an arbitrary attempt to 
force realistic circumferential stiffener geometries. The 
major effect of the additional constraint was to alter the 
dimensions of the circumferential stiffeners giving them 
more realistic values of depth and thickness. This pro- 
cedure is less direct than introducing a constraint on 
(d/t). It should be noted that the critical mode shapes 
obtained for this solution do not differ significantly 
from the critical modes shown in Table Cl. This would indi- 
cate that the 0.0063 thick by 2.0 deep ring (d/t-31 8) is 
expected to force certain critical mode shapes. This prob- 
ably is an optimistic assumption in practice. Imposition 
of the additional constraints results in a weight penalty 
of 24.9 pounds . 

To illustrate the presence of a local extremal induced 
by constraints. Case 1-1 was rerun using large initial 
constraint error weights, (4.0 x 10 6 ). The initial con- 
trol vector was taken as the final vector from the sample 
problem on page 162 of reference 9. The results obtained 
are presented in table C4. It can be seen that the solu- 
tion is constrained by the skin buckling for load case 2 
with a minimum weight of 410 pounds. This problem was 
solved above with constraint weights of 10.0 and, as noted, 
a lower minimum weight cylinder of 226 pounds was achieved. 

It should be noted that solution was simultaneously con- 
strained by four constraint functions ; gross buckling, skin 
buckling, panel buckling and longitudinal stiffener buckling. 
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The inferior solution obtained when high nominal constraint 
weights are utilized is constrained by only one constraint 
function, skin buckling, and appears to be a clear case 
of an inferior, constraint induced minima. 

Case 2-l' r Three Load Cases, Cylinder Length = 165" , 
Cylinder Radius = 60 ”. Load conditions for this problem 
are : 

1. Axial load 1400 lb/in., External pressure 

0 . 

2. Axial load 1880 lb/in.. External pressure 
-4. lb/in 2 

3. Axial load 424 lb/in.. External pressure 

.8 lb/in. 

Problem solution details are given in Table C5, Appendix C. 
Cylinder minimum weight is 387 lbs., figure 23. This com- 
pares directly to the minimum weight of 389 lbs. reported 
on page 78 of reference g. 


Case 3-1, Three Load Case, Cylinder Length = 165 " , 
Cylinder Radius = 60. " Load conditions for this problem 
are : 
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1. Axial load 2100 lb/in., External pressure 

0 . 



2. Axial load 2820 lb/in., External pressure 
-6. lb/ in 2 

3. Axial load 636 1 b/ in.. External pressure 
1.2 lb/ in 2 

Problem solution details are given in Table C-6, Appendix C. 
Cylinder minimum weight is 437 lbs., figure 23. This com- 
pares directly to the minimum weight of 445 lbs. reported 
on page 81 of reference 9. 


Case 4-0., Three Load Case, Cylinder Length = 500" , 
Cylinder Radius = 200 ". Load conditions for this problem 
are : 

1. Axial load = 2100 lb/in.. External pressure 
= 1 . lb/in 2 

2. Axial Load = 8000 lb/in., External pressure 
= -20 lb/in 2 

3. Axial load = 5000 lb/in., External pressure 
= 0 . 

This problem was solved from two nominal starting points, 
starting point 1 of reference 9 and the final solution 
obtained from starting point 1 of reference 9. Minimum 
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weights of 24252 pounds and 14332 pounds (Figure 2 3 ) > res- 
pectively, were obtained. Final cylinder details are given 
in Tables C7 and C8 of Appendix C. 

These results compare directly to a final weight of 
21300 pounds reported on page 86 of reference 9. More 
significantly, the final weight of 21300 pounds obtained 
in reference 9 is higher than the starting weight of 
16200 pounds , yet, of necessity with the Fiacco and 
McCormick procedure, both the final "minimal" weight and 
the lighter starting weight must be feasible designs. 

It should also be noted that this problem was solved 
from a second starting in reference 9. With an initial 
constraint multiplyer value of rg = 1500 a minimal weight 
of 26900 pounds was reported ;wi th an initial multiplyer 
value of rg = 375, a minimal weight of 14700 pounds was 
reported. This last solution compares favorably with the 
minimum weights of 14252 and 14332 pounds obtained by the 
AESOP exterior penalty function solution. 

Case 5-1, Three Load Case, Cylinder Length = 2000 " , 
Cylinder Radius = 200 ". Load conditions for this problem 
are : 


1 . 

Axial 

load = 

21 00 

1 b/ i n . , External 

pres sur e 


= 1.0 

1 b/ in 2 




2. 

Axial 

1 oad = 

8000 

lb/in. , External 

pressure 


= -20 

lb/in 2 
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3. Axial load = 5000 lb/in., External pressure 
= 0 

This problem was solved from a single nominal starting point 
and a minimal weight of 48097 pounds was achieved. It should 
be noted that considerable difficulty was experienced in 
obtaining this solution; the exterior penalty function con- 
straint weights were relaxed four times. It is postulated 
that the prime reason underlying these convergence diffi- 
culties was the use of relatively low constraint factors. 

In all exterior penalty function solutions of this report 
starting factors of 10. were used. The range of final 
cylinder weights covered is from 3.8 pounds to 48000 pounds. 
Clearly, the constraint penalty factors should be related 
to the payoff function values. It is suggested that studies 
using the exterior penalty function procedure should employ 
starting constraint factors approximating 1 per cent of the 
anticipated performance value. The effectiveness of alter- 
nate constraint penalty factor starting values was not in- 
vestigated in the present study. 

The one solution to Case 5 in reference 9 produced a 
minimum weight of 50000 pounds. This compares closely to 
the value of 48097 pounds above. However, in view of the 
difficulties encountered with the exterior penalty function 
solution, it is felt that a second solution to this problem 
is required to verify the optimality, or lack of optimality, 
of the result. 

Case 6-1, Single Load Case, Cylinder Length = 38" , 
Cylinder Radius = 9.55 " . The single load case is: 


Axial Load = 800 lb/in., Exterior Pressure = 0 


Two solutions were obtained corresponding to Cases 6-1 and 
6-1', pages 90 and 91 of reference 9. Minimum weights 
attained by the solutions are 3.81 and 3.70 pounds, respec- 
tively, Figure 23. This compares with 8.35 and 4.20 pounds 
in reference 9. 

Case 7-1, Single Load Case, Cylinder Length = 291" , 
Cylinder Radius = 95.5" . The single load condition was 

Axial Load = 800 lb/in.. Extremal Pressure = 0. 

This case has been discussed in detail in the section "Typical 
Optimization Algorithm Behavior in Stiffened Cylinder Design" 
where the recommended search combination produced a minimum 
cylinder weight of 682 pounds , figure 23. Solutions in 
reference 9 achieved minimum weights of 1960 pounds, 979 
pounds, and 2240 pounds, depending on the interior penalty 
function constraint factor, design variable bounds employed, 
and, in the case of the 2240 pound cylinder, the absence of 
longitudinal stiffeners. It should again be noted that the 
nominal feasible design employed (1680 pounds) was lighter 
than two of the final "mi nimum"weights reported in reference 
9. 


Case 8-1, 0, Single Load Case, Cylinder Length = 361 " , 
Cylinder Radius = 433 ". The single load condition was 

Axial Load = 12150 lb/in.. External Pressure = 0 

Starting from a nominal cylinder weighing 46840 pounds, a 
minimum weight design of 38824 pounds was attained, figure 
23. This compares to a minimum weight of 39400 pounds 
reported in reference 9. 
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Cylinder Design Summary 


It is evident from the table of results in figure 23 
that the present designs obtained by the optimization 
procedure embodied in program AESOP converge with a 
relatively high degree of reliability. The optimal 
designs obtained are all superior to thsse obtained 
previously, in some cases by a significant margin. This 
is thought to be due to three factors: 

(a) The use of multiple search algorithms 

(b) The use of an exterior penalty function 
constraint procedure 

(c) The practice of using two runs with a 
relaxation of the constraint weights 
after each run 

The solutions obtained confirm the applicability of the 
optimization techniques contained in program AESOP to 
at least one class of structural optimization problem. 
These same optimization techniques have now been applied 
successfully to problems in the fields of single vehicle 
and two -vehicle combat performance optimization, reference 
31; minimum sonic boom overpressure body shapes, reference 
32; phased array antenna design, reference 33; aerodynamic 
shaping, reference 34; liquid rocket engine combustor 
design, reference 35; and overall vehicle synthesis, 
reference 3. 


Ray Search 

The ability to perform a one-dimensional search through 
the multidimensional design space between two points P-]=a-| 
and Pg=a 2 is an integral part of current versions of the 
program of references 1 through 3. Figures 24 to 26 reveal 
the constraint and performance function behavior along such 
rays. Figure 24 illustrates function behavior between the 
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CASE 

STIFFENERS 

WEIGHT (Pounds) 

In TO e 

over 

best 

CR 

1217 
Rgsul t 

Results Obtained 
in CR-1217 

Results Obtained 
in AESOP 

Nomi nal 

Final 

Nomi na i 

Final 

1-1(1) 

Inside 

715 

231 

715 

226 

2 % 

1-1(2) 

Inside 

1000 

230 




1-1 ' 

I nside 

715 

293 




1-0 

Outside 

715 

240 




1-1,0 

Inside/ Outside 

715 

235 




1-It 

Inside* 

715 

303 




2-1 

Inside 

370 

340 




2-1 ' 

Inside 

418 

389 

418 

387 

h 

2-0 

Outside 

836 

363 



. 2-1,0 

Inside/Outside 

746 

358 




3-1 

Inside 

835 

445 

835 

437 

2 % 

3-1 ' 

I nside 

835 

490 




3-0 

Outside 

836 

468 




3-1,0 

Inside/Outside 

835 

457 




4-1 

Inside 

15900 

14600 




4-0(1) 

Outside 

16184 

21300 

16184 

1 4252 

2 % 

4-0(2) 

Outside 

44900 

26900 




4-0“ 

Outside 

44900 

1 14700 ! 



i 

4-0 (A) 

- 

- ' 

- 

21300 

1 4332 

2 % 

5-1 

Inside 

124500 

50000 

1 24500 

m 

4% 

6-1 

Inside 

13.7 

8.35 

13.7 

am 

9.5* 

6-r 

Inside 

11.8 

4.20 

11 .8 

■ 

12* 

6-0* 

Outs i de 

11.8 

4.30 


■■ 


6-1,0' 

I nside/ Outside 

11.8 

3.76 




6-Os 

Longitudinal 

12.7 

8.54 


■ 


6-Is 

Longitudinal 

12.7 

8.40 




7-1 

I n s i d e 

1680 

1960 


683 

30 % 

7-1 ' 

Inside 

1680 

979 




7-1" 

Ci rcumferenti a 1 

1680 

t 

2240 




8-1,0 

Inside/Outside 

46840 

39400 

46840 

38824 

2% 


* Reduced modulus in second load case 

** CR-1217 considers the problem to be converged 
if the value of the function is estimated to 
exceed its minimum by two per cent or less. 
Note the two per cent in half the cases. 


FIGURE 23. COMPARISON BETWEEN MINIMUM WEIGHT CYLINDER 

DESIGNS 
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4> - 48099 lbs. 

a - [.112, .1349, .01409, 2. 229, 3. 12 
(AESOP) 




FIGURE 24. CASE 5-1 RAY SEARCH 


$ » 836.18 

a f, [.03, .0335, .05009, .50166, .1 .0483, 10. 936, 1.31 1 ] 


<(> = 1682 

a = [.05, .1 , .05, 1 . 2, 8, 3] 




<(> - 836.18 * - 979 

a * [.03, .0335, .05009, .50166,1 .0483,1 0.936,1 .331 ] a = [ . 0292 , . 0441 , . 0943 , . 71 8 , . 81 , 18 . 2 ,1 .42] 
j(AESOP) (CR 1217) | 



FIGURE 26. RAY BETWEEN CR-12J7 SOLUTION AND AESOP SOLUTION, CASE 7-1 




present 48097 pound solution to Case 5-1 and the 50083 
pound solution reported in reference 9. It can be seen 
that a set of multiple-arced constraint violations lie 
between the two solutions and that an interior penalty 
function procedure would be unable to proceed along this 
multidimensional ray as a result of these violations. 
Proceeding the 50083 pound design leftwards to the 48089 
pound design constraint violations are encountered for 

(a) Rib buckling in load case 2 

(b) Skin buckling in load case 3 

(c) Skin buckling in load case 2 

(d) Rib buckling in load case 2 

(e) Gross buckling in load case 2 

Minimum weight is attained to the right of the 48089 pound 
cylinder but is accompanied by four constraint violations. 

A similar result is shown in figures 25 and 26 for 
Case 7-1. Figure 25 illustrates function behavior on the 
ray between the nominal design and the present 836 pound 
design. Figure 26 presents behavior on the ray between 
the reference 9 best 979 pound solution to this problem 
and the present 836 pound design. It can be seen from 
figure 26 that the solution shown by reference 9 is loc- 
ally constrained by skin buckling. 
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MULTIPLE EXTREMAL SEARCH PROCEDURE 


The multiple extremal search procedure proposed in ref- 
erence 1 has a straightforward basis. A multivariable 
extremal problem is solved by any of the recognized procedures, 
steepest-descent , second-order, random, or elemental pertur- 
bation techniques. With the position of an extremal known, 
a non-linear coordinate transformation which modifies the 
response surface but not the response surface topology is 
introduced. The transformation contracts an elemental "hyper 
volume" of the feasible region in the vicinity of the known 
extremal point and expands such an elemental volume at 
points far removed from the known extremal point. The feas- 
ible region boundaries are unchanged by the transformation. 

It follows directly that the "region of influence" of the 
known extremal is reduced in size and that the region of 
influence of a second extremal, if it exists, is increased 
in size. The degree of the expansion can readily be con- 
trol 1 ed . 


If the degree of the warping transformation is suffi- 
ciently high, the probability of a second multivariable 
search in the transformed coordinate space, locating the 
first extremal becomes small provided a second extremal 
point reasonably well separated from the first point exists. 
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The multiple extremal feature of program AESOP was 
demonstrated on Case 7-1 with limits on d/t (see Page 57). 
In order to use the multiple extremal feature of AESOP, 
it is first necessary to establish a warping origin in 
the parameter space. The warping origin is taken as the 
optimum design vector established in a previous run; this 
extremal is then effectively "swept out" of the response 
surface by the transformation. The exponent of the warping 
transformation in this study was taken as 2.0. The final 
control vector of the run presented in Table D8 was used 
as the warping origin. Results obtained using the warping 
transformation are presented in Table D14. (Note that this 
design was obtained with d/t limits imposed). It can be 
seen that the final weight obtained by the multiple extre- 
mal search, 838 pounds, is practically the same as that 
obtained prior to introduction of the warping transfor- 
mation, 83T pounds. 

In an additional effort to locate a relative minima, 
Case 7-1 was rerun using a 246 pound cylinder as the 
nominal. This run was made without the warping transfor- 
mation. The results of this run are presented in Table 
D15. The constraint on panel buckling is still too large 
by about three per cent; additional running would elim- 
inate this violation with no significant change in the 
cylinder weight. A final weight of 835 pounds was achieved. 
The 246 pound nominal was then rerun using the warping 
transformation centered at the final control vector of 
Table D15. The results obtained are presented in Table 
D 1 6 . 
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The weight of the lightest design obtained using 
the warping transformation is within four and one-half 
per cent of the weight of the lightest design obtained 
without application of the warping transformation. An 
additional run should significantly reduce the differ- 
ence between the two solutions. This probably indicates 
the existence of a single minimal point in the uncon- 
strained response surface. 

One objective of the present study was a demon- 
stration of the multiple extremal search on a problem 
having more than one extremal. Since this technique 
appears to be unnecessary in the case of stiffened 
cylinder design, a demonstration of the technique on a 
straightforward two-dimensional problem was undertaken. 
The problem considered is that of finding the eigen- 
values of a complex matrix. The method could readily 
be extended to the general N x N complex eigenvalue 
problem, a problem of some interest in the area of 
structural dynamics and other engineering fields. 

The 2x2 complex characteristic equation is 
( ^ 1 1 ■*" j^ii) - X a i g ^^12 
a 21 + jb 21 (a 22 + jb 22 ) - X 


= 0 

(45) 
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( 46 ) 


Here j=/^T and A is an eigenvalue of the matrix 


[A] * [a m „ + jb m J 
u J u mn mn d 

Expanding the determi nenta 1 equation (36) 

A 2 + CA + D = 0 

where 

C = A-j 1 + A 22 


(47) 


(48) 


D = 


and 


'1 1 


21 


'12 


'22 


A mn = a mn + ^ b mn 


It follows that 

A = -C*/C 2 -4L 


(49) 

(50) 

(51 ) 

This problem can readily be solved by multivariable search, 
Suppose 

A = a-j + ja 2 (52) 

The eigenvalues A are given by the points (a-j , a 2 ) which 
satisfy 

Min ( | A ( A ) j 2 ) = Min ( | A (a ] , a 2 ) | 2 ) = 0 ( 53 ) 


This problem is solved below using the matrix 


EA] = 


(1- j) 
(0+ j) 


(1 + jO) 

(0 + jO) 


(54) 
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The eigenvalues of the characteristic equation of this 
matrix are 


X-j = 1 - j 0 and X 2 = 0 - j 


(55) 


Solutions to the eigenvalue problem above were 
obtained by program AESOP using the adaptive search, an ele- 
mental perturbation technique, in conjunction with the 
"pattern" acceleration procedure. The feasibly region emp- 
loyed was defined by 


a L -j < «* 1 < a Hl 

where = aL 2 = _ 2 


a l_2 < a 2 < a ^2 

a Hi = a H2 = 2 

(56) 


The eigenvalue problem was solved from five starting points 
both with and without warping. From each point the problem 
was first solved without warping. The resulting extremal 
point was then used as the origin of a second order warping 
transformation (N = 2), and the problem was solved again from 
the same initial starting point. The second search should 
then have increased the probability of locating the second 
extremal. In all five cases the second search successfully 
found the second extremal. Figure 27 presents the results 
and the starting points. Fi gures 28a and 28b di spl ay conver- 
gence from a typical starting point. 

The example presented appears to confirm the practical- 
ity of the warping technique at least on problems of moderate 
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RESULTS OF THE FIVE MULTIPLE EXTREMAL SEARCHES 



2 

(2, 2) 

(1, 0) 

2 (Warped) 

(2, 2) 

(0,-1) 

3 

(-2, 2) 

(1, 0) 

3 (Warped) 

(-2, 2) 

(0,-1) 

4 

(2,-2) 

(0,-1) 

4 (Warped) 

(2,-2) 

_ __ __ _ _ 

(1, 0) 


5 

5 (Warped) 


(- 2 ,- 2 ) 

(- 2 ,- 2 ) 


( 0 ,- 1 ) 

( 1 , 0 ) 



















complexity (here, a fourth-order polynomial, eq. (44)). The 
regularity with which both extremals were found is surpri- 
sing. The technique merely increases the probability of 
finding a second extremal; it does not guarantee that a 
second extremal will be found if one exists. Clearly, 
further tests of the technique are in order; equally clearly, 
the multiple extremal problem is not intractable when the 
elementary technique of reference 1 can be applied so 
successful ly . 


EMPIRICAL DESIGN OF RING STIFFENED CYLINDERS 

This section contains a comparison between stiffened 
cylinder designs obtained from linear theory and ring 
stiffened cylinders obtained using empirical buckling 
data. The program used for the empirical design of ring 
stiffened shells is briefly described. In its present 
form this program does not consider longitudinally stif- 
fened shells or positive internal pressure loads. The 
comparisons are shown in figure 29 and, therefore, omit 
load conditions with internal pressure. It can be seen 
from figure 29 that the cylinders designed on an empirical 
buckling basis are between 1.50 to 3.50 times heavier 
than those designed on the basis of linear theory. These 
ratios would probably increase if the second load case of 
problems 1-1 to 5-1 were considered. 

The program uses an empirical method of analysis 
to design cylindrical or conical shell-type structures. 
Statistically determined buckling coefficients are used 
to compute the allowable loads. The program has been de- 
signed to allow buckling coefficient data to be described 
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EMPIRICAL DESIGN 


Load 

Condition 


4.34 

4.34 

4.34 

4.34 

4.34 

4.34 

13.16 

13.16 

10.53 

10.53 

6.33 

3.83 

9.5 
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6.452E-2 

.9472 

6.655E-2 

1.275E-5 

4.301E-4 


Weight 


( 1.492E-3 

703.2 

5.868E-4 ! 

489.6 

4.906E-4 

898.4 

1 . 448E-3 

629.8 

9.860E-4 1 

1038.0 

2.382E-3 

729.7 


23087.3 
29367.0 
931 23.2 
103609.9 
11,6 
1948.3 
58295.2 


LINEAR 

THEORY 


14252 


48097 

3,58 

683 

38824 


WEIGHT/W* 


FIGURE 29. COMPARISON OF CYLINDERS DESIGNED USING EMPIRICAL 
BUCKLING DATA WITH LINEAR THEORY 






















in several different ways. Most of the statistically deter- 
mined buckling data appearing in the literature can be 
described by one of the forms available in the program. For 
example, figure 30 compares the buckling coefficient for 
axial compression versus R/t used by three major aerospace 
vehicle manufacturers. The input statements of the program 
to some degree permit the user to control the method used 
to compute the allowable loads and the factor of utilization. 


Applied Loads 

Any combination of axial compression, bending, shear, 
and uniform external lateral pressure may be applied to the 
shell. Any of the above loads may be zero. 


The running load due to axial compression is computed 


by 

The maximum 
computed by 


W A = P/ [27TR. cosU)] (57) 

a 1 

value of the running load due to bending is 
W Ba - M/[irR| cos(S)] (58) 


The maximum value of the running load due to shear is 
computed by 

W Sa = V/ (ttR 1 ) + T/(2irRf) (5g) 

The running load due to uniform external lateral pressure 
is computed by 

W Qa = QR (60) 

R-j is given by 

R-j = R - L s 1 n ( £ ) / 2 c o s ( £ ) (61) 


& & 
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The Elastic Modulus at room temperature (E r j) is input to 
the program. A table is also input giving per cent of E RT 
as a function of temperature so that the elastic modulus 
at other than room temperature may be determined as a tabu 
lar function. 


Elastic Buckling Stresses 

The program has been designed to allow buckling coef 
ficient data to be described in several different ways. 

Axial Compression . The program standard option 
computes the buckling coefficient for axial compression 
as given by Seide, reference li, as 

C A = 0.606 - 0.546(1-6 1 ) (62) 

The elastic buckling stress for axial compression is then 
found to be 


for 


a A/n 


CaEcqs ( g) 



(63) 


Z < 


IT - 

^TT 


and simply supported edges 3 (see references 12, 13, and 14) 
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the buckling coefficient for axial compression is 


u - i i 36(1-P 2 )(S\Z2 

K A - 1 + —^7^ (64) 


and the elastic buckling stress for axial compression is 
then computed from 


K«TT 2 Eh 2 cos 2 (?) 

oJ n = 

rt 1 2 B 2 ( 1 - P 2 ) 


(65) 


Alternatively, in a second option, the buckling coef- 
ficient for axial compression, C ft , is input to the program. 
The user has a choice of describing the axial buckling 
coefficient as a table function or as a polynomial. The 
polynomial form adopted is 

c a ■ 2J Aa,6l! Na ' <7 (66> 

i =0 


and the elastic buckling stress for axial compression will 
be computed from 


CaEcos ( 5 ) 
°A/n = 5 


(67) 


In the third option the buckling coefficient for axial 
compression is computed from 


K A 



z > ll-.. (68) 

6C a /1 -M 2 


where C A is determined as above, and the elastic buckling 
stress for axial compression is computed as 
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V” = 


K A 7r 2 £h 2 cos 2 (5) 


(69) 


1 2B 2 ( 1 -y 2 ) 

In a final fourth program option, the buckling coef- 
ficient for axial compression is taken from reference 11 
as 


C A = 0.606 - 0.546 (1-^) + 0.16 [- h -- ° S ^ ] °‘ 3 (70) 


and the elastic buckling stress for axial compression is 
computed from 


C a E cos ( 5) 

V* " —z 


Z ^ 


6C A /l-u 2 


(71) 


For 


Z < 


6C a /Ty 2 


and simply supported edges, the buckling coefficient for 
axial compression is found from 


K 


A 


36(1-U 2 ) C A Z 2 


(72) 


and the elastic buckling stress for axial compression is 
computed from 


K fl it 2 E h 2 cos 2 (C) 

oJ n = 

A 12 B 2 ( 1 -y 2 ) 


(73) 
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Bendi ng 


The program standard option computes the buckling 
coefficient for bending from reference 11 as 

C B = 0.606 - 0.443 (1 -^) (74) 

The elastic buckling stre r for bending is then found to be 


_ C B E c °s(S) 3 ^ TOO (75) 

o R /n s 

B B Z ^ 20 

The value of o^/r\ is conservative for Z < 20. 

In a second option, the buckling coefficient for 
bending is input to the program. The user has a choice of 
describing the bending buckling coefficient as a table 
function or as a polynomial. In the polynomial form 


Nb 

C B = £ A Bi 6 1 i N b .< 7 (76) 

1-0 

and the elastic buckling stress for bending is computed as 


°B /ri 


Cb E cosU) 
3 


(77) 


The third option computes the buckling coefficient for 
bending as 

k B = 1 C B (78) 
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where Cg is determined as above, and the elastic buckling 
stress for bending is computed from 


K R t r 2 E h 2 co s 2 (5) 

o B /n = — 

B 12 B 2 ( 1 - v 2 ) 


(79) 


Shear Loads 


The standard program option computes the buckling 
coefficient for shear as given by Seide in reference 11 as 


C 


s 


0.6375 R 2 




(80) 


and the elastic buckling stress for shear is computed as 

C tt 2 E h 2 

o./n = — r. (81) 


12 B 2 (1-u 2 ) 


The value of a $ /r\ is conservative for z < 40. 


In a second option, the buckling coefficient for shear 
is input to the program. The user has a choice of describing 
the shear buckling coefficient as a table function or as a 
polynomial. In the polynomial option 


= Z 


3 / 1 


Ns 

E 

i =0 


si 


N $ .< 7 


(82) 


and the elastic buckling stress for shear is computed as 


C s tt 2 E h 2 
12 B 2 ( 1 - u 2 ) 


(83) 
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Uniform External Lateral Pressure 

In the standard program option, the buckling coeffi- 
cient for uniform external lateral pressure is taken from 
reference 11 and the elastic buckling stress for uniform 
external lateral pressure is computed as 


K p tt 2 E h 2 
12 B 2 (1-y 2 ) 


(84) 


In a second option the buckling coefficient for 
uniform external lateral pressure is input to the program. 
The user has a choice of describing the pressure buckling 
coefficient as a table function or as a polynomial. In 
the polynomial option 


N p 

K P ■ E v z1 ; V ' 1 < 85 > 

i = 0 


and the elastic buckling stress for uniform external 
lateral pressure is computed as 


K p tt 2 E h 2 
12 B 2 (1-u 2 ) 


( 86 ) 
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Plasticity Correction 


In the standard program option, the elastic buckling 
stresses are not corrected for the effects of plasticity. 
The allowable buckling stresses are set equal to the elas- 
tic buckling stresses 


a = a/n (87) 

In a second program option the elastic buckling stresses 
are corrected for the effects of plasticity. The following 
procedure is used to compute the allowable buckling stresses. 


If 

a/n ^ a uE 

then a = a uA 

(88) 

If 

a/n ^ a L 

then a = a/n 

(89) 

If 

a L < a/n < a uE 

then a will be 

computed 


a = (T-T.) + a. (90) 

2 ~ 1 1 1 

where 

T i T ^ T ^ 

and 

N 

a i = ^ A1 ^o/n) 1 N < 4 
i =0 

N 

a _ ^ A 2.j (a/n)^ N 4 4 
i =0 



Ratio of Applied to Allowable Loads 


The allowable buckling stresses described above are 
used to compute the allowable loads and the ratio of 
applied to allowable load as follows: 

The allowable axial compressive loads are given by 


w. 


h a, 


(91) 


and 


Pqr = 2tt w R- j cos(^) 


(92) 


The ratio of applied axial load to allowable axial load is 

P (93) 


Ri 


CR 


The allowable bending loads are given by 


W B ^ a B 


(94) 


and 


'CR 


= TT Wg R-J 2 cos(£) 


(95) 


The ratio of applied bending load to allowable axial load 
i s 


R B ” M 


M 


CR 


(96) 


The allowable shear loads are given by 


w s = h a $ 


(97) 
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The ratio of applied shear load to the allowable shear load 
i s 


R 


s 


"s» /w s 


(98) 


The allowable uniform external lateral pressure load is 
given by 


w 


P 


= h Op 


(99) 


and 


g CR • Wp/R (100) 

The ratio of applied pressure load to allowable pressure 
load is 


R 


P 


- <WcR 


( 101 ) 


Interaction of Combined Loads 


In the standard program option the exponents used in the 
interaction equation are given by 


Y 


1 + 0.7(|) - 0.04 (f) 




1 0 


y = 2.0 


f* 10 


Y 


P 


= T 


= Y B = Y s 


1.0 


( 102 ) 



Alternately, the exponents y^, y B> y g , y, and y p may 
be input as data. 


The factor of utilization is computed from 


1.0 = 


FU Y B Y S 

<tt> * (ir> + <§*■> 


,Rd, Y P 


(S 2 ) 


(103) 


or optionally as 


1.0 = |(^) YA + + (^-) YS + (^) YP 004) 


Ys ,Rn Y P 


The margin of safety is found from 

M.S. = {(1/U) - 1} 


005) 


Frame Stiffness Requirements 

The required frame stiffness, Ely, is given by 
Timoshenko, reference 15, for cylinders subject to uniform 
external lateral pressure and axial compressive force as: 
( NOTE : The spacing of the frames is assumed to be small 

compared to the radius of the shell) 


a 


1 


ly H=^l 
B h R 2 


Where Iy is the effective Iy 
when the shell is sub- 
jected to the design 
loads (106) 


a-j is obtained from 

0-j + C^ot + = + 0g<j>2 


(107) 
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where 


C ] = sX 4 

C 2 = X 6 (X z +2n z ) + sX 2 n 2 [2U 2 -l ) 2 +2(n 2 -l ) 2 + 5X 2 n 2 - 2] 

C 3 = ( n 2 - 1 ) 2 [X 4 +s ( 2X 2 +n 2 ) n 2 ] 

= X 4 n 2 + s(2X 2 + n 2 )n 4 - s ( 3 X 2 + n 2 )n 2 

C 5 = X 6 + sX 2 n 2 (2X 2 + n 2 + 1) (108) 

Set m = 1 and let n = 2, 3, 4, . . . , so that the maximum 
value of a-j may be determined. 

£n the above equati ons the parameters a, X, s hy, , <j> 2 , and 
P were determined from 


h 


X = 


mitR 


S - hy mp- 


hy = h + (Ay/B) 


4 - QRO-k ) 

*1 Eli 


and 


P (1-y 2 ) 
2tt R Eh 


* 

P 


2M 

R 


+ P 


(109) 


For the value of computed above to be the valid, the 
following must be true: 


U = 1.0 


(no) 
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Hess, reference 16, gives a graphical method for 
solving the above equations for . In the present study 
the equations are solved numerically. 


Wei ghts 

Shell weight is given by 

WT S = 2ttRL hp (111) 

For frame weight, assume "z" section frames and flange 
length to thickness ratio of 10 and a frame web flat 
height of 20t. The frame area is approximately 

A F = 40 4^74 012) 

Normally, the program uses the criteria that the two end 
frames must provide "fixed" edges assuming that an end 
ring stiffness to intermediate frame stiffness ratio of 
20 will provide "fixed" ends. With this assumption, the 
total weight of the frames is approximately given by 



The end frame weights are omitted in the present study. 
Total weight is given by 

WT T0TAL = WT F + WT s (114) 


(113) 
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CONCLUSION 


It is apparent that modern multivariable optimization 
techniques are capable of solving the class of structural 
design problems considered in this report. Optimal designs 
are achieved routinely in a reasonably efficient manner; 
each solution obtained required the expenditure of approx- 
imately three minutes' time on a high speed large-scale 
digital computer (the CDC 6600) using a single load case. 

The study results lend support to the view that topo- 
graphically complex constrained optimization problems can 
be more reliably solved by sequential application of 
several search algorithms than by the repeated application 
of a single search algorithm. This premise underlies the 
optimizing program AESOP which was employed in the study. 

It is considered significant that the search combination 
utilized throughout the major portion of the study consisted 
of a random technique, an elementary perturbation technique, 
and two straightforward search acceleration techniques. 

Over the spectrum of problems considered in this study, no 
advantage was perceived in the use of more organized tech- 
niques such as gradient methods or second-order methods. 

In point of fact, experience in the solution of both the 
present problems and a variety of optimization problems in 
other fields tends to support the view that as response 
surface complexity increases, the selected search procedure 
should be weighted towards the use of techniques such as 
the random ray procedure. 
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It would appear that the exterior penalty function 
approach is well suited to solution of multivariable 
optimization problems which involve non-convex constraint 
boundaries. This characteristic of the exterior penalty 
function technique is dependent on an ability to penetrate 
constraint boundaries in pursuit of Improved performance. 
With non-convex constraints the improved performance can 
often be retained provided the constraint boundary subse- 
quently enters the region of improved performance. 

Successful exploitation of the exterior penalty function 
technique in the presence of non-convex constraint boun- 
daries is somewhat dependent on the use of adaptively 
determined constraint weighting factors and a willingness 
to restart the solution using relaxed weights following 
convergence to an initial constrained optimum. 

No evidence of multimodal behavior in the unconstrained 
response function (cylinder weight) surface itself was 
detected. When the unconstrained response surface itself 
possesses more than one extremal, the search techniques 
applied here can be combined with a true multiple extremal 
search procedure such as the topographically invariant 
warping of program AESOP. 

It is well known that a linear buckling analysis will 
result in an unconservative design. This point is clearly 
demonstrated by the present study; for minimum weight 
designs are obtained by the application of both linear and 
empirical buckling criteria. Cylinder weights obtained by 
the two approaches differ by factors as high as 3.5 even 
though in some cases the critical load case was not con- 
sidered in the empirically based design. Inclusion of this 
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load case can only lead to higher weight empirically 
designed cylinders. It is accordingly recommended 
that the input of practical buckling criteria be con- 
sidered in future studies of the present type. It is 
also recommended that future studies incorporate 
realistic geometries on stiffener members. 

The successful application of multivariable 
search techniques to the stiffened cylinder design 
problem encourages their further application to struc- 
tural design. Ultimately, one seeks a method capable 
of practical application to large-scale general purpose 
structural analysis programs such as the National 
Aeronautics and Space Administration NASTRAN program. 
Further development of the multivariable search pro- 
cedures may be required before such an approach becomes 
practical with today's computers. On the next generation 
of computers such as the CDC 7600 or possibly the STAR 
computers present techniques would appear capable of 
optimal structural definition through general purpose 
codes provided efficient multiple analysis techniques 
are employed. 
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APPENDIX A 


Development of the Analysis of the Stiffened Cylinder 
A.l Introduction 

In this appendix all the equations needed to analyze the stiffened 
cylinder are presented. These include the overall buckling analysis of 
the cylinder as well as the buckling, stress and yield analyses of the skin 
and stiffeners. 

It is well known that there is a large discrepancy between the 
buckling failure loads for monocoque cylinders which are predicted by 
classical buckling theory and the failure loads obtained in tests. However, 
it has been found recently that this is not necessarily the case for 
stiffened cylinders, reference 17. Linear theory is used here but it has 
been found that this may not apply in some cases, reference 18. The 
importance of including the effect of eccentricity of the stiffeners 
has been pointed out both experimentally, reference 19, and analytically, 
references 20, 21, and 22. Earlier investigators have also treated this 
effect analytically, references 23 and 24. In the analysis used here, 
eccentricity effects are included. This analysis follows closely that 
of Flugge in reference 24. 


A. 2 Stress-Strain Relations 

The skin of the cylinder is assumed to be in a biaxial state of stress. 
The axes of elastic symmetry are in the longitudinal and circumferential 
directions. The x axis is in the longitudinal direction and the <f> axis is 
in the circumferential direction. With these assumptions the stress-strain 
relations in the sheet are 
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The stiffeners are assumed to be in a uniaxial state of stress so 
that the stress-strain relations are 
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in the longitudinal and circumferential stiffeners respectively. 


A. 3 Strain-Displacement Relations 

The reference surface of the shell is taken as the midsurface of the 
skin. With the z axis taken positive inward from the reference surface 
and u, v, and w being the displacements of the reference surface respectively 
in the positive x, <j>, and z coordinate directions, the strain displacement 
relations are taken to be 
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where ex» e^,, and y x< f, are the strains at a point in the shell e x and 
e<j> are assumed to be continuous in the skin and x and <j> stiffeners, 
respectively. These relations may be derived in a geometric manner as 
done by Fliigge, pg. 212 of reference 24, or by reducing the linear 
three-dimensional strain displacement relations in cylindrical coordi- 
nates, reference 25. The latter is done by assuming the displacements 
vary linearly with the depth of the shell, reference 26, and by setting 
the transverse shear strains and the extensional strain in the z direction 
to zero. 


The displacements of a point in the cylinder corresponding to these 
strain midsurface displacements are 
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(see Figure Al). 

The rotations of the normal used in the above displacements are 
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The relative rotations per unit length are then 
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A. 4 Force Resultants 

The force and moment resultants per unit length are obtained by 
performing the appropriate integrations of the stresses over the thickness 
of the skin and then adding to these the corresponding force and moment 
resultants per unit length in the stiffeners. The force and moment 

resultants per unit length in the stiffeners are obtained by dividing the 
resultant forces and moments by the stiffener spacings. 

The extensional forces and bending moments in the stiffeners are 
obtained by performing the appropriate integrations of the extensional 
stresses in the stiffeners over the areas of the stiffeners. The stiffeners 
are assumed to carry no shear load; so they have no contribution to the shear 
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force resultants, but they are assumed to have a twisting moment resultant. 
The angle of twist is assumed to be the same as that of the normal to the 
skin. The torsional stiffnesses of the stiffeners are obtained from an 
approximate curve for data given by Crandal and Dahl, reference 27, for 
torsion of bars of rectangular cross section. Thus, the force resultants 
are obtained by substituting the strain displacement relations (A3) into 
the stress strain relations (Al) and (A2) , then substituting the resulting 
stress displacement relations into the following formulas and performing 
the integrations: 
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The above expressions apply for stiffeners on the inside of the cylinder. 

For stiffeners on the outside of the cylinder the limits of integration on 

the stiffener integrals, the second terms, must be changed to go from 

- (d + t/2) to - t/2 and - (d^ + t/2) to - t/2 (see Figures 17 and A2) 

e x and are the angles of twist of the normal to the skin, given in 

section A. 3. and J, are the section constants for a rectangular cross- 
X 9 

section in torsion. These correspond to a polar moment of inertia and are 
approximate by the expression, 

J = c ab^ b <_ a 

where c is given by 

c = - 0.285 e "°- 49(a/b) + 0.316 
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and a and b are the cross sectional dimensions of the stiffener, t x and d x , 

and t A and d,; b is taken as the dimension of smaller magnitude. After 

<p v 

making the substitutions described above, performing the integrations, and 
neglecting terms of the order of the thickness of the skin divided by the 
radius and square of the depth of the stiffeners divided by the square of 
the radius with respect to 1, the force and moment resultants can be written: 
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where the constants in the above expressions are given in terms of the 
material properties and the dimensions of the stiffened cylinder by 
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The effects of the eccentricity of the stiffeners are seen in the 
3 

terms e x , e^, and , which have a positive sign when the stiffeners 
are on the inside and a negative sign when the stiffeners are on the outside. 


A. 5 Pre~buckle Forces and Stresses 

It is assumed that when the cylinder is loaded there is a uniform 
change in length and a uniform change in radius. This implies that u, v, 
and w are independent of <p; w and v are independent of x; and that u is 
a linear function of x. Applying these assumptions to the force displace' 
ment relations (A. 7), the forces in the cylinder are 
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By substituting these into the equilibrium equations, reference 24, 
page 209, the internal forces may be obtained in terms of the applied loads. 
The result is 
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where N is the applied axial compression load per unit length of circum- 
ference, and p is the applied external pressure per unit surface area. 

With the assumptions about the prebuckled deformation, the midsurface 
prebuckle strains are obtained from the strain displacement relations as 
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By equating the expressions for the force resultants in terms of the dis- 
placements with the values of the force resultants in terms of the external 
forces and identifying the strains, the following expressions for the mid- 
surface strains are obtained, after neglecting terms of the order of the 
depth’ of the stiffener divided by the radius with respect to one: 
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Substituting these into the stress-strain relations (A.l) 
the expressions for the stresses in the skin are obtained: 
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The expressions for the stresses in the ribs neglecting terms involving 
the depth of the stiffener divided by the radius with respect to one are 
obtained by multiplying the stiffener modulus by the corresponding value 
of the midsurface prebuckle strain. These are 
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A. 6 Buckling of the Cylinder and Skin 

An expression for the critical buckling load of the cylinder is 
obtained in terms of two integer parameters representing the buckling mode 
shape. The lowest buckling load is then obtained by searching the buckling 
loads obtained from a large number of possible mode shapes. 

The expression for the buckling load is obtained from the determinant 
of a set of homogeneous equations. These are obtained by substituting into 
the buckling equilibrium equations, in terms of displacements, an assumed 
solution which satisfies these equations and simple support boundary condi- 
tions. The displacement functions contain the two integer parameters 
representing the mode shape and arbitrary constants. 

The buckling equilibrium equations are obtained in terms of displace- 
ments, by substituting the force resultants, in terms of the displacements, 
into the buckling equilibrium equations in terms of forces. The buckling 
equilibrium equations used are those given by Flugge, reference 24, page 422, 
but contain only the buckling force terms recommended by Hedgepeth and Hall, 
reference 17, page 9. With the changes required because of the different 
coordinate system used here these equations are 
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where N is the applied axial compressive force per unit length and p is the 
applied external radial pressure per unit area. 

After substituting the force displacement relations (A7) into these 
equations, the buckling equilibrium equations in terms of displacement are 
obtained. These can be written in the form: 
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The assumed displacements which satisfy the above equations and the simple 
support boundary conditions are 


u = A sin n<t> cos xx 

v = B cos n<(> sin xx 

w = C sin n4> sin xx 


where for the complete cylinder 

X = f, m *1,2,... 

n=n n - 0,1,2,... 

and L is the length of the cylinder. 
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For a cylindrical plate (the skin between stiffeners) 
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After the displacements are substituted into the displacement buckling 
equilibrium equations, these equations can be written in the form: 
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where the C's are given by 
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Since the values of the applied loads are known, a ratio between the 
axial load and pressure can be calculated. Letting 


p = a N 

and setting to zero the determinant of the coefficients of A, B, and C in 
the last set of equations the expression for the critical axial load is 
obtained. This is 


, N v = - F + V^B 2 - 4 A C ' 

R s2 cr zK 

where 

A = C u (R x 4 + x 2 (n 2 - 1) a) R 

B = [(C-ji C 22 - C 12 ) + (C-ji C 33 - C-| 3 C 31 )] R X 

+ (C n C 22 - C ]2 2 ) a (n 2 - 1) 
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c 


C 11 C 22 C 33 


+ c 12 c 23 c 31 + c 13 c 12 c 23 


" C 12 2 C 33 ‘ C 13 C 3l C 22 “ C 23 2 C 11 


For each combination of the parameters m and n there are two possible values 
of (N/H s2 ) cr . The one which has to be used as critical is the one with 
smallest magnitude and having the same sign as the applied load N. The 
critical buckling load is obtained by finding (N/H s2 ) cr for a large number 
of values of both m and n and then selecting the lowest magnitude value 
out of all of these. 

For the special case when N = 0 the critical pressure must be found. 
This is given by 


(£— ) 

H s2 cr 


- C 

(('ll C 22 " C 12 2 ) ^ 2 " ^ 
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The above analysis is used for gross buckling, panel buckling, and 
sheet buckling. For gross buckling all the constants are calculate as 
given in the C‘s and the full length of the cylinder is used. For panel 
buckling, the terms which contain the properties of the circumferential 
stiffeners are set to zero and the length of the cylinder is taken as the 
circumferential ring spacing. For skin buckling, all terms containing 
stiffener properties are set to zero, n is changed to apply to a cylin- 
drical plate with a width of the longitudinal stiffener spacing and the 
length of the cylinder is again taken as the length between circumferential 
sti ffeners. 


A. 7 Longitudinal Stiffener Buckling 

The critical buckling stress for the longitudinal stiffeners is 
obtained by applying a solution for the critical buckling stress of a flat 
rectangular plate to several different possible assumed modes of buckling 
of the stiffener. In all the possible assumed modes the longitudinal 
stiffener is assumed to be simply supported on three edges and free on the 
fourth. The critical buckling stress for such a flat plate is given by 
Bleich, reference 28, as 


A15 


(A23) 


cr 


c 


12(l-v ) d 


, 2 

[(f) + 


0.425] 


The notation has been changed here; t is the thickness of tne plate (i.e. 
the width of the stiffener), d the width of the plate (i.e. the depth of 
portion of the stiffener under consideration), and a the length of 
stiffener under consideration. 

The first failure mode to which this expression is applied is in the 
situation where the circumferential stiffeners are either on the opposite 
side of the cylinder from the longitudinal ones or where they are non- 
existent. In this case d is taken as d x , the full depth of the stiffener, 
and sl is taken as L, the full length of the cylinder. 

The second mode is where the circumferential stiffeners are on the 
same side of the cylinder as the longitudinal ones and are the deepest. 

In this case the critical buckling stress of the longitudinals is taken as 
that of a plate with depth d x , the full depth of the longitudinal stiffen- 
ers, and a length $, x , the length between circumferential stiffeners. 

The third mode is where the circumferential stiffeners are on the 
same side as the longitudinal ones but are not as deep. In this case one 
would expect the stiffener to buckle in a manner coupling the material 
between the circumferential stiffeners with the material above the circum- 
ferential stiffeners. To obtain an estimation of the critical buckling 
stress two cases are considered. One assumes that the portion of the 
material between the circumferential stiffeners does not buckle but the 
outstanding portion does. In this case the formula is applied to a plate 
of the dimensions of the depth of the outstanding portion and the length of 
the entire cylinder (d x - d^ by L). The other case assumes that the 
material between the circumferential stiffeners does buckle with the out- 
standing portion of the longitudinal stiffener but that the circumferential 
stiffeners force nodes in the buckling of the longitudinal and these notes 
occur at the location of the circumferential stiffeners. The buckling 
stress in this case is taken as that for a plate and d equal to the full 
depth, d x , of the stiffener and i equal to i, x , the circumferential 
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stiffener spacing. This is the same as the case where the circumferential 
stiffeners are deeper than the longitudinals. 

A. 8 Circumferential Stiffener Buckling 

Similar situations are encountered with the buckling of the circum- 
ferential stiffeners as with the buckling of the longitudinal stiffeners. 
Here, however, an additional mode of buckling Is encountered (see TablesAl 
and A2). The external stiffeners not only can buckle when they are 
compressed, but due to their curvature can also buckle when they are 
expanded. An expression for the critical circumferential strain in the 
skin of the cylinder, or at the edge of the stiffener, is obtained (in 
Section A. 9) by doing an assumed mode solution of the buckling problem. 

This expression is 


<(icr 


t . 2 

( 4 > 


(— ■ (i + u ± > x 

12(l-v ) 1+2? + C 


(A24) 


2 2 

(1 + 2n 2 (l-v)) c + (2n 2 (2-v)-l) + (n 2 - 1) ? 3 / 3 + ... 

[ „ , ] 


i * (2^H) 1 7 (2^, £ 2 7 


where d is the depth of the stiffener portion in question; and c is the 
ratio of the stiffener depth, d, to the radius of the unsupported edge of 
the stiffener; c is a positive number if the stiffener is inside and nega- 
tive if the stiffener is outside; and n is the number of full waves in the 
circumferential direction. 

With the circumferential stiffeners on the inside of the cylinder, c 
positive, the value of e^ cr is negative for all values of n and increases 
in magnitude as n increases. This means that inside circumferential 
stiffeners can buckle only when the cylinder contracts under load. 

With the circumferential stiffeners inside the cylinder and the 
longitudinals outside or non-existent the critical buckling value for e 

<pcr 

is obtained with n = 0. 
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With both the circumferential and the longitudinal stiffeners inside 
the cylinder and with the longitudinal stiffeners deeper than the circum- 
ferential ones, the circumferential stiffeners are physically restrained 
from buckling into a smaller number of half waves than the number of spaces 
between longitudinal stiffeners. Since e^ cr increases in magnitude as n 
increases, the critical buckling value for this situation is obtained by 
using for n the number of spaces in half the circumference of the cylinder. 

In the situation with the circumferential stiffeners deeper than the 
longitudinal ones two values of e^ cr are obtained. One is for the un- 
supported portion of the circumferential with d 5 d • d„ and n = 0. The 

<p X 

other is obtained as above, for the supported stiffeners, for the full depth 
of the stiffen er assuming that nodes are forced at the locations of the 
longitudinal stiffeners. This is similar to the case of the longitudinal 
stiffeners. 

With the circumferential stiffeners outside, t, negative, e^ cr is 

positive for small values of n and increases in magnitude as n increases. 

When n becomes large enough e,^„ becomes negative and then as n increases 

<*>cr 

the magnitude of e^ cr decreases while the value remains negative. The 
magnitude of E^ cr decreases until for some value of n a minimum is obtained. 
Thus the circumferential stiffener can buckle for small values of n when the 
cylinder expands, e positive, or can buckle for large values of n when 

Cf) w i 

the cylinder contracts, e (J>cr negative. 

For the case of external circumferential stiffeners with the longi- 
tudinal stiffeners inside or non-existent, the critical positive value of 

e . is obtained with n = 0, and the critical negative value is found by 

per 

searching for the lowest magnitude negative value of e ^ cr . 

With the longitudinal stiffeners also on the outside and deeper than 
the circumferential the circumferential stiffeners are again physically 
restrained from buckling into a smaller number of half waves than the 
number of spaces between the longitudinal stiffeners. A value of e is 

(pC ( 

calculated for n equal to the number of spaces in half the ci rcumference 
n = ttR/a^ . If this value is positive then this is the critical value for 

an expansion of the cylinder. If it is negative there is no critical value 

for an expansion of the cylinder. Several possibilities exist for the 
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negative buckling value. If the above value of e^ cr is positive then the 
negative value is the one given by the minimum magnitude value found for the 
unsupported stiffener, since this value has a larger number of circumfer- 
ential waves than spaces between stiffeners. If the value for n = irR /£ 
is negative then there is a choice between this value and the value 
minimum |e A I. The one which has the larger value of n is used. The 

9C1 

reasons for this are as follows: if n = is the largest then a 

smaller n is physically impossible; if the n for minimum [e | is larger 
then this gives smallest |e^ cr | for e^ cr negative and is physically 
possible. For the case of all external stiffeners with the circumferential 
ones having the greater depth the problem is again split into two parts, one 
an unsupported circumferential stiffener with depth d^ - d x and the other 
a stiffener with the full depth d^ , assuming nodes at the location of the 
longitudinal stiffeners. Values are then obtained for each case in a 
manner similar to that described above for external stiffeners. Two 
results are then obtained and compared to find the critical value. 

In this treatment n is considered as a continuous variable instead 
of integer as it actually is and no arguments about the compatibility of the 
mode shapes are made. Introducing these restrictions would increase the 
buckling values so that the treatment used is conservative. 

This buckling solution does not apply where the ci rcumferential 
stiffener is thick compared with its depth. This is because the assumed 
simply supported boundary condition does not apply. In situations where 
the outstanding portion of the stiffener has a depth to thickness ratio of 
less than ten the yield limit is substituted for the buckling limit. 

A. 9 Solution of the Circumferential Stiffener Critical Buckling 
Strain 

An approximate solution is obtained for the buckling of a circular 
plate with a large hole in it. The plate is assumed to be simply supported 
at one edge and free at the other (see Figure A3). The simply supported 
edge is the edge which attaches to the cylinder and thus must have the same 
displacements as the cylinder. The critical buckling parameter is taken as 
the tangential strain on the simply supported edge. The solution is 
obtained using an assumed mode variational method. 



The variational formulation of the problem of elastic stability is 
given by Novozhilov, reference 26, page 173, as 

<5[ A^ ] = <5 


In the cases in which the initial stress state can be determined using 
classical theor 
coordinates as 


( 2 ) 

classical theory, this is such a case, A v ' is given in cylindrical 
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The primes denote the buckle state and zero the initial state. By the same 
type of procedure as used for the derivation of the strain displacement 
relations in the cylinder the above strain displacement relations can be 
reduced to strain midsurface displacement relations. Thus, by assuming 
the displacements vary linearly with depth, the transverse shear strains are 
zero, and the normal strain is zero, the strain displacement relations reduce 
to 

2 

, _ 3U 1 _ 3 W 1 

£r " ar 2 

2 

i _ 1_ av 1 u 1 z A a w 1 . aw 1 \ 

e 6 r 39 r ” r 'r ..2 ar ' 

o y 

2 

1 = v 1 , 1 au 1 _ Zz_ i a w 1 _ 1_ aw 1 \ 

e re ar ~ r r as r 'ar ae " r ae ' 

i 1 aw 1 , aw 1 , av 1 , v_L _ 1 au 1 
“r r ae ’ “e ar * “z ar r ~ r ae 


Note that the displacements in these expressions- are the midsurface displace- 
ments and not the displacements of a point as in the previous expressions, 
also u 1 , v 1 , and w 1 are the displacements in the r, e, and z directions and 
are not the same as the u, v, and w for the cylinder. 


The stresses in the plate prior to buckling are given by Timoshenko, 
reference 25, page 59, in terms of a radially inward pressure at the outer 


edge as 
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P is the radial pressure. These may be transformed to be in terms of the 

tangential strain at the edge r = b, F„ , by using the stress strain 

u 

relations to solve for P in terms of e 0 and substituting the result 
into the above expressions for the stresses 


b2 ‘a - 7» % 
r a 2 + b 2 + v(a 2 - b 2 ) 


( A25) 


h 2 e B * 7> 

6 a 2 + b 2 + v(a 2 - b 2 ) 


The following set of buckling displacements satisfy the displacement 
boundary conditions: 


u‘ = v' = 0 


w‘ = A (b-r) sin ne 


These displacements are then substituted into the strain-displacement rela- 
tions and the resulting expressions along with the prebuckle stresses are 

(21 

then substituted into the expression for A' . When the integration is 
carried out and s[A^] is set to zero (6 r| 2 ^ is zero for the problem 
since the forces are constant on one edge and the displacements are 
constant on the other), reference 26, pg. 172, the following expression is 
obtained for the critical £ e : 

- _ t 2 r a 2 + b 2 + v(a 2 - b 2 )-, 

*• ' ' mu7) : F -3 

(„ 2 -l> 2 1„| + 2„ 2 (l-„ 2 ) „ 2 

[ 2 * 3 (A26) 

l-2n /k 2 2^ /- u 2^2 , ^ 2 f 2 -|^\ b 

— 7£— (b - a ) + (b n + a (n -1 ) ) In — 
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Now is set equal to ? and the expression is expanded in terms of 

cl ^ 

this quantity. The result of doing this and setting the critical e Q equal 
to the critical strain in the cylinder is the expression: 


<*cr 
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Specialization of this solution in two limiting cases, for which solutions 
are available in the literature, is given in Appendix B. 


A. 10 Yield Failure 

The principal stresses in the skin are given by a vn and a (A13). 

xp 9P 

It is assumed that the yield criterion for the cylindrical shell skin 
material is of the following form, reference 29. 
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where 
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the longitudinal tension yield stress in the skin if a X p > 0 
the longitudinal compression yield stress in the skin if 
a xp < 0 

the circumferential tension yield stress in the skin if 


the circumferential compression yield stress in the skin 

if cr < 0 

<PP 


constant defining yield envelope in first quadrant a > 0 

xp 

and « > 0 

4>P 
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^ae = K CT constant defining yield envelope in second quadrant a xp < 0 

and 0 > 0 

4>P 

= K CC constant defining yield envelope in third quadrant a xp < 0 

and a. n < 0 

< PP 

K aS = k TC constant defining yield envelope in fourth quadrant o xp > 0 

and o.„ < 0 

<t>P 

For the case of an isotropic material that behaves identically in tension 
and compression with yield stress oq^ Eqs. A27 when specialized by the 
following substitutions: 


a xOT 


< 


a3 


a xOC = a ^0T = a <pOC °0D 
1 


reduce to the distortion energy yield criterion 


2 , 2 2 
a xp ‘ a xP%p a <j,p - a 0D 
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The stiffeners are in a uniaxial state of stress so the stresses (A14) 
must satisfy the yield conditions: 


a xS0C - °xsp - a xS0T 
a ,#,S0C - %sp - °<#>S0T 


(A29) 


where the subscript x refers to the longitudinal stiffener; <j> refers to 
the circumferential stiffener; 0 refers to yield; C refers to compression; 
and T refers to tension. 
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TABLE AT. SELECTION OF CIRCUMFERENTIAL STIFFENER BUCKLING MODE (ej - CONTRACTION 
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FIGURE A2. FORCE RESULTANTS 
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APPENDIX B 


Verification of the Circumferential Stiffener 
Buckling Solution 


Two limiting cases for e r are obtained. One is for the washer 
mode of buckling with no circumferential waves. The other case is for a 
large number of waves. These cases are checked with existing solutions. 

With n = 0 and neglecting terms involving the depth of the stiffener 
divided by the radius, c, with respect to one, the expression for the washer 
mode is 


£ <pcr 
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— *r (2)t 
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With n large and again neglecting c with respect to one the express 
ion for e , „„ is 

( p Cl 
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Substituting the expression for the critical strain in the washer 
mode (Bl) into the expression for the stress in the radial direction (A25) 
and obtaining the value at the supported boundary, r = b, the following 
expression is obtained: 


-E._ t 2 

°rcr = ^ 2 $ « (2 C ) 

rCr 12 (1 - v 2 ) d 

Writing c as d/R and making the definition D = E^ s t 2 /1 2(1 - v 2 ) this 
expression is 


The negative sign signifies that the stress is compressive. The above value 
agrees closely with the exact solution for the axisymmetric case done by 


Bl 


Meissner, reference 30, who obtain a value for the coefficient 1.86 instead 
of 2. 

By substituting the value of e ecr for n large (B2) into the expres- 
sion for the tangential stress, a Q , (A25) at the supported edge, r = b, the 
following expression is obtained: 

~ E ,c + 2 2 2 

°e = ^ — o~ (J) 0(1 - v) + S-JL _ ) ( 2 ) 

6 12(1 - v 2 ) d 2 

n is expressed in terms of the half wavelength, i, as ttR/a . Substituting 
for e and n the following expression is obtained: 


ecr 


:hk . J (t) 2 (( 1) 2 + 

12 (1 - V 2 ) d * ir 2 


If the value \> = .3 is used this becomes 


^ecr 



+ 2 H 2 
(?) «f> 


0.425) 


which is the same as the expression for the critical buckling stress for a 
rectangular plate (see Section A. 7). 
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APPENDIX C - DETAILED SOLUTIONS FROM EXTERIOR PENALTY FUNCTION METHOD 

Case I -I 




fc x 

% 

d x 

d 4 

*x 

% 

W 

Final 

.02179 

.03047 

.00629 

. 3608 

2.0 

7.2379 

.8768 

225.93 

Initial 

.099 

.06 

.06 

.5 

.5 

6.0 

3.0 

714.92 

U. B. 

.5 

.5 

.5 

2.0 

2.0 

o 

o 

■ — i 

10.0 


L. B. 

.0001 

.0001 

.0001 

.0001 

.0001 

.oooi 

.0001 



f/f 


cr 


L.C. 

1 

2 

3 

G. B. 

.7893 

1.0044 

.2450 

P.B. 

.8545 

.9918 

.3782 

S . B. 

.9153 

1.0026 

.3380 

LRB. 

.7195 

1.0046 

.2102 

CRBU 

0 

0 

0 

CKBL 

-8427. 

-17874. 

-1240. 

S.Y. 

.4041 

.5863 

.1164 

LRYU 

-.4057 

-.5664 

-.1185 

LRYL 

.4057 

.5664 

.1185 

CRYU 

.1261 

. 2676 

.1857 

CRYL 

-.1261 

-.2676 

-.1857 

LOADS 

LC 

N 

P 


1 

700. 

0 . 


2 

940, 

-2, 


3 

212. 

.4 



L « 165. R = 60. 


,v 

Aluminum = .101 

E = 10 x 10 6 v - .333 


a y =* 50,000. 


Wave Numbers 


LC 

Gross 

Panel 

Skin 

M 

14 

1 

8 

1 

N 

9 

31 

1 

M 

o 

14 

1 

10 

4 

N 

9 

24 

1 

M 

13 

1 

7 

j 

N 

9 

97 
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TABLE Cl 






C*se 1 - I (B) 


o 

rv> 




fc x 

% 

d x 

d * 


s 

w 

Final 

. 01982, 2_ 

.03097 

.00851 

.36516 

1,91928 

7.1545 

.7971 

227.72 

Initial 

.02000 

.03120 1 

.00907 

.36500 

1.92000 

7.1500 

.7920 

230.39 

U. B. 

.5 

.5 

.5 

2.0 

2.0 1 

10.0 

10.0 


L. B. 

.000001 

.000001 

.000001 

.000001 

.000001 : 

. 000001 
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cr 


L.C. 

1 

2 

3 

G.B. 

,776.8. 

. 9999 

.2449 

P.B. 

.7724 

.9084 

.3665 

S.B. 

. 9247 

. 9999 

.3481 

LRB. 

.7170 

.9999 

.2097 J 

CRBU 

0 

0 

0 

CRBL 

-4189 . 

-9150. 

-563.8 

S.Y. 

. 4058 

.5891 

.1172 

LRYU 

-.4080 

-.5690 

-.1194 

LRYL 

. 4080 

. 5690 

.1194 

CRYU 

.1232 

. 2691 

; .1658 

CRYL 

-.1232 

-.2691 

-.1658 
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LC 

N 

P 


1 

700. 

0 


2 
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-2. 


3 

212 

.4 



L » 165. 


R = 60. 


Aluminum 
E = 10 x 10 6 


Y = .101 

v = .333 0 =, 50,000. 

y 


Wave Numbers 


LC 

Gross 

Panel 

Skin 

M 
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N 

13 

1 

9 

9 

30 

1 

M 

o 

14 

1 

10 

4 

N 

9 

23 

1 

M 

o 

13 

i 

7 

J 

N 

9 

115 

1 


TABLE C2 






Case 1 - I (A) 

(with |crbl| < l.o) 



c s 

fc x 

% 

d x 
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Z x 

% 

W 

Final 

.02345 

.03198 

.04912 

.3936 

. 8847 
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.5 

.5 

o 

VD 
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o 
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.5 

.5 

.5 
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O 

O 
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L. B. 

.0001 

. 0001 

.0001 

. 0001 

. 0001 

.0001 

.0001 



o 

Cl) 


f/f 


cr 


L.C. 

i 

2 

3 

G.B. 

.8022 

1.0032 

.2894 

P.B. 

.8868 

1.0020 

. 3993 

S.B. 

.8979 

1.0039 

.3258 

LRB. 

. 7206 

1.0020 

.2114 

CRBU 

0 

0 

0 

CRBL 

-.1428 

-.3024 

-.02114 

S.Y. 

.3727 

. 5338 

.1083 

LRYU 

-.3761 

-.5229 

-.1103 

LRYL 

.3761 

. 5229 

.1103 

CRYU 

.1045 

.2212 

; .01546 

CRYL 

-.1045 

-.2212 

-.01546 

LOADS 

LC 

N 

P 


1 

700. 

0. 


2 

940. 

-2. 


3 

212. 

.4 



L = 165. R = 60. 


Aluminum 

E = 10 X 10 6 V 


.101 

.333 


Oy - 50,000. 


Wave Numbers 


LC 

Gross 

Panel 

Skin 

M 

11 

1 

8 

1 

N 

10 

29 

1 

M 

o 

12 

1 

10 

z 

N 

10 

22 

1 

M 

o 

1 

1 

7 

J 

N 

5 

84 

1 


TABLE C3 






o 

4 ^ 


CASE I 

Using Large Constraint Weights 4 x lo 6 



fc x 

% 

d x 

d » 

04649 

.1052 | 

RSM 

i .3341 

. 3225 

05455 


.1381 

.3238 

.3659 

5 

.5 

.5 

2.0 

2.0 

019 

.05 

.05 

0.0 

0.0 


10.0 


.05 


2.48 


2.57 


10.0 


.05 


410.15 


489.78 


8833 

.9537 

.8130 

4155 

1 .5065 

.1423 

8725 

. 9994 

.3144 

0291 

.0407 

.0085 


E - 10x10” 


V = .333 


Wave Numbers 
LC I Gross 1 Panel 


CRBL 

-.0698 

-.1395 

-.0120 

S.Y. 

.2278 

. 3271 

.0658 

LRYU 

-.2289 

-.3193 

-.0669 

LRYL 

.2289 



.3193 

.0120 

CRYU 

.0698 ! 

.1395 ; 

.0120 

CRYL 

-.0698 

-.1395 

-.0120 


LC 

N 

P 

1 

700 

0 

2 

940 

-2 

3 

l 212 

0.4 


TABLE C4 

















































iSE 2-1' 


dx 

. a* . 

*X 

& ^ 

W 



m.nn. 

1.2078 

387.40 

.519 

.129 

9.68 

1.16 

417.57 

2.0 

2.0 

10.0 

10.0 


■2.0 

iLi 

.05 

.05 



L = 165. R = 60. 

y = .101 

E = 10.0 x 10 6 v = .333 cr y = 50,000. 


WAVE NUMBERS 


LC 

GROSS 

PANEL 

SKIN 

M 

mm 

1 

8 

1 




N 

9 

24 

1 

i 

10 

1 

9 

B9 

9 

19 

1 

M 

1 

1 

6 

3 

N 

4 

68 

1 


TABLE C5 

































































































































Case 4-0 Starting Point 1 



fc s 

fc x 

% 

d x 

d * 

*x 

-e- 

w 

Final 

.1438 

.1820 

.2158 

-1.750 

-2.3028 

39.96 

4.6015 

14252.6 

Initial 

.15 

.20 

.489 

BBH 

-1.89 

35.0 

4.58 

16184.3 

U. B. 

1. 

20. 

40. 

mam 

4. 

40. 

20. 


L. B. 

.00001 


. 00001 

-4. 

-4. 

.00001 

.00001 



o 


L = 500 


R s 200 


L.C. 

1 

2 

3 

G.B. 

.9904 

.9720 

.8009 

P.B. 

.5593 

. 9841 

. 8732 

S.B. 

.3265 

.5621 

.6643 

LRB. 

.2224 

1.008 

.5457 

CRBU 

.0358 

.7257 

.1443 

CRBL 

-.0335 

-.6792 

-.1350 

S.Y. 

.1869 

1.036 

. 4645 

LRYU 

-.1901 

-.8618 

-.4665 

LRYL 

.1901 

.8618 

.4665 

CRYU 

. 3584 

.7257 

.1443 

CRYL 

-.3584 

-.7257 

-.1443 


LOADS 


LC 

N 

p 

1 

2100. 

1. 

2 

8000. 

o 

o 

CNJ 

3 

5000. 

0 


•Y .101 

Aluminum 

E = 10 x 10 6 V — .333 » 50,000. 


Wave Numbers 


LC 

Gross 

Panel 

Skin 

M 

1 

1 

7 

1 

N 

6 

48 

1 

M 

o 

13 

1 

14 

N 

24 

6 

1 

M 

3 

1 

8 

J 

N 

9 

27 

1 


TABLE C7 





























Starting Point B 



TABLE C8 

















































CASE 5-1 



fc s 

*x 

S 

d x 

d » 



W 

Final 

.11224 

.13488 

.01409 

2.2286 

3.1201 

1.4105 

6.4853 

48097 

Initial 

.25 

.25 

.3 

2.0 

10. 

25. 

4. 

124500 

U. Be 

.5 

20. 

. . _.! 

1 . 

5. 

10.1 

100. 

10. 


L. B. 

.1 

.02 

.01 

0.0 

0.0 

1. 

1 . 



L.C. 

1 

2 

3 

G.B. 

1.0028 

1.0028 

0.8713 

P.B. 

P . 0050 

0.0190 

0.0119 

S.B. 

3.2018 

0.6969 

0.4732 

LRB 

3,1220 

0.5407 

0.2981 

CRBU 

D.O 

0.0 

0.0 

CRBL 

-95.738 

-18282.9 

-3743.0 

S.Y. 

3.1730 

0.8872 

0.4259 

LRYU 

-0.1762 

-0.7805 

-0.4304 

LRYL 

3.1762 

0.7805 

0.4304 

CRYU 

3.0294 

0.5614 

0.1149 

CRYL 

-0.0294 

-0.5614 

-0.1149 


LOADS 


LC 

N 

r 

P 

1 

2100. 

1 . 

2 

8000. 

-20. 

3 

5000. 

0 . 


L ■ 2000. R m 200. 
a y = 72000. Y * .101 
E = 10 . 5xl0 6 V = .333 


Have Numbers 


LC 

Gross 

Panel 


M 

1 

1 

3 

Bl 

N 

3 

0 

1 

M 

41 

3 

1 

2 




N 

8 

0 

1 

M 

31 

3 

1 

J 




N 

9 

0 

1 


TABLE C9 

















































CIO 


CASE 6-1 



t s 



dx 

d<j) 

*x 

& (j) 

— 

Final 

Initial 

U.B. 

L.B. 

.013634 

.015188 

.000533 

.29760 

.24387 

.19178 

2.000 

3.8135 

.04 

.04 

.04 

.25 

.25 

1.0 

1.0 

13.723 

1.0 

1.0 

1.0 

2.0 

2.0 

5.0 

2.0 


.0000001 

.0000001 

.0000001 

2.0 

zZJl 

.0000001 

.0000001 


L = 38.0 R= 9.55 


L.C. 

1 

G.B. 

1.00660 

P.B. 

.01564 

S.B. 

.59987 

LRB 

.69918 

CRBIJ 

0 . 

CRBL 

-1012.93 

S.Y. 

.99959 

LRYU 

-1.00215 

LRYL 

1.00215 

crytj 

.31958 

CRYL 

-.31958 

LOADS 


L.C. 

N 

P 

1 > — 1 

1 

800. 

0 . 


Y= .101 

E = 10.5x107 y= .333 o y = 50000. 

WAVE NUMBERS 


L.C. 

GROSS 

PANEL 

SKIN 

M 

6 

4 

1 

1 




N 

8 

0 

1 


TABLE CIO 















































CASE 7-1 



fc s 

fc x 

s 

d x 

d * 

% 

X 

Him 

■*_ 

Final 


.0276 


.3879 

20. 0__ 



Buncm 1 

Initial 

.05 

.i 

.05 

1.0 

2.0 

8.0 

3.0 

1681.74 | 

U. B. 

10.0 

ZSHM 

10.0 

20.0 

20.0 

50.0 

20.0 


L. B. 

.0000001 

.0000001 

.0000001 

.00000001 

.0000001 

.0000001 

.0000001 



L. C . 

1 

2 

3 

G.B. 

1.0028 



_ _ __ 

P.B. 

.2173 

— 

— 

S.B. 

1.0051 




LRB. 

1.0071 

— 

— 

CRBU 

0.0 

— 

— 

CRBL 

-1.5xl0 10 

— 

— 

S. Y. 

.4145 



— 

ini 

-.4146 

— 


LRYL 

.4146 

— 

— 

CRYU 

.1375 

— 



CRYL 

0.1375 

— 

— 


LOADS 


LC 

N 

P 

1 

800 

0.0 

2 

i 

i 

i 

i 

i 

— 

3 1 

— 

— 


L = 291 R = 95.5 

Y „ .101 

E = 10 x 10 6 V = .333 


■ 50000 


Wave Numbers 


LC 

Gross 

Panel 

Skin 

M 

1 

N 

27 

i 

2 

6 

62 

1 

M 

2 

N 


„ 

„ 




M 

3 

N 



__ 


— 

— 

— 


TABLE Cl 2 















































I-C 


1 


G.B. 
P.B. 
S • B. 
LRB . 
CRBU 
CRBL 
S~Y. 
LRYU 
LRYL 
CRYU 
CRYL 


1,0041 

.4992 

1.0028 

1.0057 

0.0 

-3.57x10 
. 6384 
- . 6408 
.6408 
. 2000 
-.2000 


3 


co 


LOADS 
LC N 

JL 12150 

2 


3 


CASE 8-1 ,0 



TABLE Cl 3 
















































APPENDIX D 

SOLUTION TO CASE 7-1 USING SEARCH COMBINATION (9, 2, 9, 2,5, 2, 3) 




Sc 

% 

d x 


l 

X 

% 

w 

Final 

.03044 

.0276 

.000022 

.3879 

20.0 

3.229 

1.3162 

682.54 

Initial 

0.05 

0.1 

0.05 

1.0 

2.0 

8.0 

3 * 0 

1681.7 

U. B. 

10.0 

10.0 

10.0 

20.0 

20.0 

50.0 

.20.0 


L. B. 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 



L.C. 

1 

G.B. 

1 . 002*8 

P.B. 

.2173 

S.B. 

1.Q051 

LRB 

1.0071 

CRBU 

0 

CRBL 

-1.5E+10 

S.Y. 

.4145 

LRYU 

41 4'6 — 

LRYL 

.4146 

CRYU 

.1375 

CRYL 

-.1375 


NUMOPT « 
METHOP - 


LOADS 


LC 

N 

P 

1 

800 

0 



L - 

291.0 

R “ 95.5 

at. 

tf 

910000.0 

7 

Oy * 

50000. 

r - o.ioi 

« 

tX' 


9,2, 9, 2, 5, 2, 3 

E * 

10.5E6 

V = .333 

14.1 


Wave Numbers 


LC 

Gross 

Panel 

Skin 

M 

1 

27 

1 

2 

N 





6 

62 

1 


SEARCH TECHNIQUES USED: 

Random Ray (Method 9) 
Pattern (Method 2) 

Creeping (Method 5) 
Magnification (Method 3) 


TABLE Dli 





SOLUTION TO CASE 7-1 USING SEARCH COMBINATION (5, 2, 5, 2) 


O 

ro 



fc s 


s 

d x 

d * 


% 

w 

Pinal 

.10904 

.00959 

.00020 

.000008 

2.00737 

2.46866 

20.0 

1925.84 

Initial 

0.05 

0.1 

0.05 

1.0 

2.0 

8.0 

3.0 

1681.7 

U. B. 

10.0 

10.0 

10.0 

20.0 

20.0 

50.0 

20.0 


,L. B. 

l.OE-7 i 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 



L.C. 

— I T 

G.B. 

1.00367 

P.B. 

.37317 

S.B. 

.36318 

LRB 

1.1E-9 

CRBU 

0 

CRBL 

6.1E+8 

S.Y. 

.1467 

LRYU 

-.1467 

LRYL 

.1467 

CRYU 

.04879 

CRYL 

-.04879 


NUMOPT - * 

METHOP - 5 , 2 , 5 . 2 


LOADS 


LC 

N 

P 

1 

800 

0 


L - 

291.0 

R « 

95.5 

d<f> 

ty * looooeo 

°y “ 

50000. 

Y “ 

0.101 

dx 

£7 ■ .00089 

E « 

10.5E6 

v = 

0.333 


Wave Numbers 


LC 

Gross 

Panel 

Skin 

M 

1 

N 

48 

1 

1 

14 

0 

i 

1 


SEARCH TECHNIQUES USED: 

Creeping (Method 5) 
Pattern (Method 2) 


TABLE D2 

















SOLUTION TO CASE 7-1 USING SEARCH COMBINATION (6,2, 6, 2) 



S 


% 

d x 

d * 


% 

w 

Final 

.05785 

1 . OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

12.4752 

1020.25 

Initial 

0.05 

0.1 

0.05 

1.0 

2.0 

8.0 

3.0 

1681.7 

U. B. 

10.0 

10.0 

10.0 

20.0 

20.0 

50.0 

20.0 


.L. B. 

1 . OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

1 


O 

CO 


L.C. 

1 

G.B. 

3.616-2* 

P. B. 

.00369 

S.B. 

.00369 

LRB 

.00100 

CRBU 

0 

CRBL 

-. 0921 

S. Y. 

. 2765 

LRYU 

-.2765 

LRYL 

.0921 

CRYU 

-.0921 

CRYL 



LOADS 


NUMOPT « 4 
METHOP - 6 ' 2 ' 6 ' 2 


* Excessive constraint 
violation 


LC 

N 

p 

1 

800 

0 


L « 291.0 
a y - 50000. 
E « 10.5E6 


LC 


Panel 

Skin 

M 

1 

N 

15 

11 

11 

I 

30 

ipper Limit 
100 

15 


R « 95.5 -4a “ 1 

Y « 0.101 

Ay 

V = 0.333 || - 1 

Wave Numbers 


SEARCH TECHNIQUES USED: 

Quadratic (Method 6 ) 

Pattern (Method 2) 

NOTE: The results presented are- those 

obtained from the first run. The 
second run was aborted due to a 
singular matrix. 












SOLUTION TO CASE 7-1 USING SEARCH COMBINATION (7, 2, 7, 2) 


■ 

% 


% 

d x 

d * 

*x 


w 


.02260 

.03022 

.00794 

.41874 

2.54229 

7.45131 

.93042 

684.44 


0.05 

0,1 

0.05 

1.0 

2.0 

8.0 

3.0 

1681.7 

U. B. 

10.0 

10.0 

10.0 

20.0 

20.0 

50.0 

20.0 


L. B. 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 



L*C« 

1 

G«B. 

i.04 as * 

P. B. 

.6944 

S .B. 

1.0178 

LRB 

1.0123 

CRBU 

0 

CRBL 

-1.09E+6 

S . Y. 

.4362 

LRYU 

- . 4386 

LRYL 

.4386 

CRYU 

.1320 

CRYL 

-.1329 


LOADS 


NUMOPT - H 
METHOP - 7 ' 2 ' 7 ' 2 


Excessive constraint 
violation 


LC 

N 

P 

1 

800 

0 


L * 291.0 
Qy * 50000. 
E - 10.5E6 


d4> 

R - 95.5 * 321.0 

Y - 0.101 

dx 

V = 0.333 tx « 13.8 


Wave Numbers 


LC 

Gross 

Panel 

Skin 

M 

1 

N 

19 

1 

8 

10 

40 

1 


SEARCH TECHNIQUES USED: 

Davidon (Method 7) 
Pattern (Method 2) 


TABLE D4 




















SOLUTION TO CASE 7-1 USING SEARCH COMBINATION (4, 2, 4, 2) 





% 

d * 


% 

X 

s 

W 

Final 

.10321 

l.OE-7 

.00001 

l.OE-7 

7.32029 

4.3547 

2.89986 

1820.4 

Initial 

0.05 

0.1 

0.05 

1.0 

2.0 

8.0 

3.0 

1681.7 

U. B. 

10.0 

10.0 

10.0 

20.0 

20.0 

50.0 

20.0 


L. B. | 

l.OE-7 

l.OE-7 ! 

l.OE-7 

l.OE-7 

l.OE-7' 

l.OE-7 

l.OE-7 



o 

cn 


L.C. 

1 

G.B. 

1.117.4* 

P.B. 

1.0110 

S.B. 

.1467 

LRB 

.0019 

CRBU 

0 

CRBL 

-8.9E+11 

S.Y. 

.1550 

LRYU 

-.1550 

LRYL 

.1550 

CRYU 

.0516 

CRYL 

-.0516 


LOADS 


NUMOPT 

METHOP 


4 

4 t 2 r 4 t 2 


291.0 R - 95.5 

50000 Y * 0.101 

10.5E6 M = 0.333 

Wave Nimbers 


•d$ 

t<r 


• 732029.0 


fix _ , 
HT" 1 


LC 

Gross 

Panel 

Skin 

M 

1 

Jpper Limi 
50 

1 

2 

N 





12 

0 

1 


LC 

N 

P 

i 

800 

0 


Excessive Constraint Violation 

TABLE D5. 


SEARCH TECHNIQUES USED: 

Steepest-Descerlt (Method 4 ) 
Pattern (Method 2) 























a 

CTi 


SOLUTION TO CASE 7-1 USING. SEARCH COMBINATION (5,2,3) 




fc x 

% 

d x 

** 

*x 

% 

w 

Final 

.07795 

.13221 

.03189 

.01691 

.67382 

2.71116 

.86511 

1558.30 

Initial 

0.05 

0.1 

0.05 

1.0 

2.0 

8.0 

3.0 

1681.7 

U. B. 

10.0 

10.0 

10.0 

20.0 1 

20.0 

i 

50.0 

20.0 


L. B. 

1.0E-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 



L.C. 

1 

G.B. 

1.7832 * 

P.B. 

1.0227 

S.B. 

. 03301 

LRB 

. 00003 

CRBU 

0 

CRBL 

“94.37 

S'. Y* 

.1957 

LRYU 

-.1967 

LRYL 

.1967 

CRYU 

. UbUU / 

CRYL 

-.06007 


NUMOPT « 2 
METHOP » 5 


L - 

291.0 

R * 95.5 

t* 

21.1 

Oy « 

50000. 

y - o.ioi 


E - 

10.5E6 

V = 0.333 

dx m 
tx 

.0128 


Wave Numbers 


LC 

Gross 

Panel 

Skin 

M 

1 

N 

Jpper Iini-t 
50 

1 

3 

15 

0 

1 


SEARCH TECHNIQUES USED: 

Creeping (Method 5) 
Pattern (Method 2) 


* Excessive constraint Violation 


LOADS 


LC 

N 

P 

i 

800 

0 


TABLE D6 
























SOLUTION TO CASE 7-1 USING SEARCH COMBINATION (1,2) 




Sc 

' s 

d 

X 


l 

X 

% 

w 

Final 

. 00941 

. 08609 

l.OE-7 

.96933 

20.0 

.03294 

16.679 

255.08 

Initial 

0.05 

0.1 

0.05 

1.0 

2.0 

8.0 

3.0 

1681.7 

U. B. 

10.0 

10.0 

10.0 

20.0 

20.0 

50.0 

20.0 


-L. B. 

T 1.0E-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

i 


L.C. 

i 

G.B. 

1.130&* 

P.B. 

.0227 

S.B. 

.01480 

LRB 

.00084 

CRBU 

0 

CRBL 

-1.8E+17 

S.Y. | 

1.1095 

LRYU 

-1.110* 

LRYL 

1.110* 

CRYU 

.3675 

CRYL 

-.3675 


LOADS 


NUMOPT « * 
METHOP - 1,2 


♦Excessive constraint 
violation 


LC 

N 

P 

i 

800 

0 


L ■ 

291.0 

R - 95.5 

||* ■ 20.OE+7 

a v * 

50000. 

Y - 0.101 

jf 

E - 

10.5E6 

V = 0 .333 

s- 


Wave Numbers 


LC 

Gross 

Panel 

Skin 

M 

1 

N 

13 

6 

14 

7 

19 

3 


SEARCH TECHNIQUES USED: 

Sectioning (Method 1 ) 
Pattern (Method 2) 


o 

^4 


TABLE D7 






SOLUTION TO CASE 7-1 USING SEARCH COMBINATION (9 , 2 , 9 , 2 , 5 , 2 , 3) 
o (Limits imposed on d/t) 



fc s 

fc x 

% 

** 

d * 

Z 

X 

% 

w 

Final 

0.030 

0.0335 

0.05009 

0.50166 

1.0483 

10.936 

1.3110 

831.21 

Initial 

0.05 

0.1 

0.05 

1.0 

2.0 

8.0 

3.0 

1681.7 

U. B. 

10.0 

10.0 

10.0 | 

20.0 

20.0 

50.0 

O 

o 

jjHHBf | 

L. B. 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 



L.C. 

1 

G.B. 

1. 0046* 

P.B. 

1.0058 

S.B. 

1.0072 

LRB 

1.0028 

CRBU 

0 

CRBL 

-57.114 

s.y. 

0.3695 

LRYU 

-0.3721 

LRYL 

0.3721 

CRYU 

0.1083 

CRYL 

-0.1083 


NUMOPT * 
METHOP - 


LOADS 


LC 

N 

P 

1 

800 

0 


L 

7 

9,2 9, 2, 5, 2, 3 E 


291.0 

50000. 

10.5E6 


R - 95.5 

Y » 0.101 

V = 0.333 


d4 

t<j> 


21.0 



15.0 


Wave Numbers 


LC 

Gross 

Panel 

Skin 

M 

r 

15 

1 

8 

N 





13 

35 

1 


SEARCH TECHNIQUES USED: 

Random Ray (Method 9) 
Pattern (Method 2) 
Creeping (Method 5) 
Magnification (Method 3) 


TABLE D8 























SOLUTION TO CASE 7-1 USING SEARCH COMBINATION (5, 2, 5, 2) 
(Limits imposed on d/t) 




Sc 

% 

d x 

d * 

X 

S 

w 


.02446 

.03211 

.06172 

.46184 

1.2766 

10.280 

1.0132 

819.04 


0.05 

0.1 

0.05 

1.0 

2.0 

8.0 

3.0 

1681.7 

0. B. 

10.0 

10.0 

10.0 

20.0 

20.0 

50.0 

20.0 


L. B. 

l.OE-7 

L ... . 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 



L.C. 

1 

G.B. 

1.0044. 

P.B. 

1.0050 

S.B. 

1.0050 

LRB 

1.0033 

CRBU 

0 

CRBL 

-54.339 

s.y. 

.3979 

LRYU 

-.4024 

LRYL 

. 4024 

CRYU 

.1048 

CRYL 

-.1048 


NUMOPT « 4 
METHOP - 5,2, 5,2 


L - 291.0 
<7y * 50000. 

E « 10.5E6 


R « 95.5 20.7 

Y - 0.101 

v = 0.333 ^. = 14.4 


Wave Numbers 


LC 

Gross 

Panel 

Skin 

M 

1 

N 

16 

1 

9 

11 

35 

1 

1 


SEARCH TECHNIQUES USED: 

Creeping (Method 5) 
Pattern (Method 2) 


LOADS 


LC 

K 

P 

i 

800 

0 


Table D9. 


o 























DIO 


SOLUTION TO CASE 7-1 USING SEARCH COMBINATION (6, 2, 6, 2) 
(Limits inposed on d/t) 




*x 

% 

d x 

d * 

*x 

S 

W 

Final 

.02544 

.02891 

.05625 

.44930 

1.3713 

9 . 5067 

1.0375 

806.85 

Initial 

0.05 

0.1 

0.05 

1.0 

2.0 

8.0 

3.0 

1681.7 

U. B. 

10.0 i 

10.0 

10.0 I 

20.0 

20.0 

50.0 

20.0 


.L. B. 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 



L.C. 

1 

G.B. 

L . 0036 

P.B. 

L . 0056 

S.B. 

L . 0025 

LRB 

1.2045* 

CRBU 

0 

CRBL 

-75.909 

s.y. 

.40922 

LRYU 

-.4139 

LRYL 

.4139 

CRYU 

.1074 

CRYL 

-.1074 


NUMOPT « 4 
METHOP - 6, 2, 6, 2 


* Excessive Constraint 
LOADS Violation 


L - 291.0 
Cy - 50000. 
E « 10.5E6 


R m 95.5 * 24.4 

Y m 0 • 1 0 1 
■ dx 

v = 0.333 37 - 15.6 


Wave Numbers 


LC 

Gross 

Panel 

Skin 

M 

1 

N 

18 

1 

8 

10 

37 

1 


SEARCH TECHNIQUES USED: 

Quadratic (Method 6) 
Pattern (Method 2) 


LC 

N 

P 

i 

800 

0 


Table DIO. 






SOLUTION TO CASE 7-1 USING SEARCH COMBINATION (7,2) 
(Limits imposed on d/t) 




Sc 

% 

d x 

** 

*x 


w 

Final 

.03718 

.03222 

.03948 

.50308 

.90802 

8.3231 

1.7334 

894.04 

Initial 

0.05 

0.1 

0.05 

1.0 

2.0 

8.0 

o 

m 

1681.7 

U. B. 

10.0 

10.0 

o 

o 

20.0 

20.0 

50.0 

20.0 


-L. B. 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 



L.C. 

1 

G.B. 

1.0072 

P.B. 

.6931 

S.B. 

1.0044 

LRB 

.9977 

CRBU 

0 

CRBL 

-64.618 

S.Y. 

. 3389 

LRYU 

-.3407 

LRYL 

.3407 

CRYU 

.1029 

CRYL 

-.1029 


LOADS 


NUMOPT * 
METHOP « 


LC 

N 

P 

1 

800 

0 


2 

7,2 


L 

°y 

E 


291.0 

R * 95,5 

■d$ 

t# 

- 20.8 

50000. 

Y * 0.101 


10.5E6 

V = 0.333 

d* 

tx 

- 14.1 


Wave Numbers 


LC 

Gross 

Panel 

Skin 

M 

1 

N 

16 

1 

5 

14 

40 

1 


SEARCH TECHNIQUES USED: 
Davidon (Method 7 ) 
Pattern (Method 2) 


Table Dll. 





















D12 


SOLUTION TO CASE 7-1 USING SEARCH COMBINATION (4,2) 
(Limits imposed on d/t) 




*x 

s 

*x 

d * 

A 

X 

% 

w 

Final 

.02970 

.03158 ' 

.04092 

.46631 

1.1744 

8.3364 

1.2821 

824.52 

Initial 

0.05 

0.1 

in 

o 

o 

1,0 

2.0 

8.0 

3.0 

1681.7 

U. B. 

io.o I 

10.0 

10.0 

20.0 

20.0 I 

50.0 

20.0 


La B. 

l.OE-7 i 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 



L.C. 

1 

G.B. 

.99903 

P.B. 

.73365 

S.B. 

1.0029 

LRB 

1.0053 

CRBU 

0 

CRBL 

-111.53 

S.Y. 

. 3803 

LRYU 

-. 3834 

LRYL 

.3834 

CRYU 

.1089 

CRYL 

-.1089 


NUMOPT - 2 
METHOP - 4,2 


LOADS 


LC 

N 

P 

1 

800 

0 


L - 

291.0 

R« 95.5 

t <P 

28.3 

Oy * 

50000. 

Y - 0.101 



E - 

10.5E6 

V = 0.333 

tx 

14.75 


Wave Numbers 


LC 

Gross 

Panel 

Skin 

M 

1 

N 

17 

1 

6 

12 i 

39 

1 


SEARCH TECHNIQUES USED: 

Steepest-Descent (Method 4) 
Pattern (Method 2) 


Table D12 





























SOLUTION TO CASE 7-1 USING SEARCH COMBINATION (1,2) 
(Limits imposed on d/t) 




fc x 

% 

d x 

d * 

*x 

K JL 
$ 

w 

Final 

.02561 

.03075 

.05682 

. 43385 

1.1801 

9.3976 

.97670 

813.46 

Initial 

0.05 

0.1 

0.05 

1.0 

2.0 

8.0 

3.0 

1681.7 

U. B. 

10.0 

10.0 

10.0 

20.0 

20.0 

50.0 

20.0 


L. B. 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 



L.C. 

1 

G.B. 

1.069a* 

P.Bp 

.9811 

S.B. 

. 8421 

LRB 

.9622 

CRBU 

0 

CRBL 

-59.967 

s.y. 

^ . 3969 

LRYU 

-.4010 

LRYL 

.4010 

CRYU 

.1070 

CRYL 

-.1070 


LOADS 


NUMOPT *2 
METHOP *1,2 


^Excessive constraint 
violation 


LC 

N 

P 

i 

800 

0 


L * 

291.0 

R ■ 95.5 

a 

20.8 

Gy * 

50,000. 

Y * 0.101 


E * 

10.5E6 

V « 0.333 

fe- 

14.4 


Wave Numbers 


LC 

Gross 

Panel 

Skin 

M 

1 

17 

1 

9 

N 



■ ■ 


11 

37 

PPM 


SEARCH TECHNIQUES USED: 

Sectioning (Method 1) 
Pattern (Method 2) 


Table D13 





















SOLUTION TO CASE 7-1 USING AESOP WARPING TRANSFORMATION 


o 

-P* 





S 

d 

X 


l 

X 

i 

$ 

W 

Final 

.02899 

.03203 

.04921 

.47911 

.98 

6.9709 

1.1208 

837.76 

Initial . 

0.05 

0.1 

1.05 

1.0 

2.0 

8.0 

3.0 

1681.7 

U. B. 

10.0 

10.0 

10.0 

20.0 

20.0 

50.0 

20.0 


L. B. 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 


WARPAL 

0.03 

0.0335 

0.05009 

0.50166 

1.0483 

10.936 

1.3110 



L.C. 

1 

G.B. 

. 9951 

P.B. 

. 9644 

S.B. 

. 9795 

LRB 

.9986 

CRBU 

0 

CRBL 

-56.66 

S.Y. 

. 3669 

LRYU 

-.3695 

LRYL 

. 3695 

CRYU 

.1075 

CRYL 

-.1075 


LOADS 


NUMOPT «= 
METHOP = 


LC 

N 

P 


800 

0 


7 


9, 2, 9,2, 5, 2,3 


= 291.0 

R - 95.5 

% * 19 - 8 

t<t> 

= 50000. 

Y = 0.101 

= 10.5E6 

V = 0.333 

^ =*15.0 
tx 


Wave Numbers 


LC 

Gross 

Panel 

Skin 

M 

1 

N 

15 

1 

8 

13 

36 

1 


Search Techniques Used: 

Pattern (Method 2) 

Random Ray (Method 9) 
Creeping (Method 5) 
Magnification (Method 3) 


TABLE D14 



















D15 


SOLUTION TO CASE 7—1 USING SEARCH COMBINATION (9 /2 / 9 / 2 /5 / 2 / 3) 



t 3 


% 


% 


% 

• I 

Final 

.02979 

.03376 

.05185 

.50544 

1.0381 

11.239 

1.3178 


Initial 

.01 

.01 

.01 

.2 

.2 

1.0 

1.0 

246.74 | 

U. B. 

10.0 

10.0 

10.0 

20.0 

20.0 

50.0 

20.0 


L. B. 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 



L.C. 

1 

G.B. 

1.0018 

P.B. 

1.0329 

S.B. 

1.0061 

LRB 

1.0006 

CRBU 

0 

CRBL 

-51.736 

S.Y. j 

.3677 

LRYU 

-.3703 

LRYL 

. 3703 

CRYU 

.10789 

CRYL 

-.10789 


NUMOPT ■ 
METHOP - 


LOADS 


LC 

N 

P 

1 

800 

0 


L - 

7 °Y “ 

9,2,9 r 2 ( 5,2,3 E " 


291.0 R * 

50000. f ■ 

10.5E6 v = 


95.5 

d* 

20.0 

0.101 

dx 

15.0 

0.333 

tx - 


Wave Numbers 


LC 

Gross 

Panel 

Skin 

M 

1 

N 

15 

1 

8 

13 

35 

1 


SEARCH TECHNIQUES USED: 
Random Ray (Method 9) 
Pattern (Method 2) 
Creeping (Method 5) 
Magnification (Method 3) 


TABLE D15 



































SOLUTION TO CASE 7-T USING AESOP WARPING TRANSFORMATION 


a 

O') 



fc s 

5c 

s 

d x 

d * 

l 

X 

l 

* 

W 

Final 


. 10350 


1.3206 

.41088 


20.0 

867.48 

Initial 

.01 

.01 

.01 

.2 

.2 

1.0 

1.0 

246.74 

U. B. 

10.0 

10.0 

10.0 

20.0 

20.0 

50.0 

20.0 


L. B. 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 

l.OE-7 



L.C. 

1 

G. B. 

1.0045 

P. B. 

.0079 

S . B. 

1.0029 

LRB 

.4275 

CRBU 

0.0 

CRBL 

-1.304 

S.Y. 

.4500 

LRYU 

-.4568 

LRYL 

.4568 

CRYU 

.1011 

CRYL 

-.1011 


291.0 

R = 

95.5 

50,000 

Y = 

.101 

10.5E6 

v 5=5 

.333 


Wave Numbers 


I LC 


Panel 

Skin 

■ 

6 

2 

1 

11 

0 

1 


LOADS 


LC 

N 

P 

1 

800 

0 


Table D16 








































APPENDIX E 


PROGRAM STRUCTURE AND DATA INPUT/OUTPUT DESCRIPTION 


Figure El defines the overlay structure for the 
program. The central memory core storage requirement 
is 70000 q locations. All data required by the program 
is input through two namelist data blocks. The first 
data block, n CR1217 n , is for the analysis data; the 
second data block, IAESOP, is for the optimization data. 

For each complete analysis the following output is 
obtained : 

CLT These are the critical buckling load 

values divided by H 2 for the modes saved 
for the skin for the last load condition. 

SMS The values of m saved, starting with all 

the values for gross buckling for all 
load conditions , followed by panel buck- 
ling for all load conditions, followed 
by skin buckling for all load conditions . 

SNS Same as SMS but for the values of n which 

are saved. 

CRITICAL LOADS Each line contains the critical buckling 
load for gross, panel, and skin buckling 
for one load condition; successive load 
conditions are on successive lines . 

MODE SHAPES Same order as above giving the values of 

m and n. 

LRS Stress in the longitudinal rib for each 

load condition . 

CRS Stress in the circumferential rib for each 

load condition. 

DES Actual value of distortion energy stress 

squared for each load condition. 


El 


EBU 


Critical strain value, circumferential rib, 
for an expansion of the cylinder, for each 
load condition, 

EBL Critical strain value, circumferential rib 

for a contraction of the cylinder. 

LRCB Critical buckling stress for the longitudi- 

nal rib, for each load condition. 

BEU Logical variables signifying the existence 

of a critical strain EBU, T for True, F for 
false. 

BEL Same as above for EBL. 

BLR Same as above for LRCB. 

EPA Actual value of circumferential strain for 

each load condition. 

TS, TX, TY, DX, DY, LX, LY These correspond to t s , t x , 

t(j) # d x , ^(j) f A<j) • 

AX, AY Areas of the longitudinal and circumferen- 

tial stiffeners , respectively. 


The following eleven lines of output are the ratios 


of actual 
values , in 

values of the behavior variables to the critical 
columns for each load condition. 

G.B. 

Gross buckling 

P.B. 

Panel buckling 

S.B. 

Skin buckling 

LRB 

Longitudinal stiffener buckling 

CRBU 

Circumferential stiffener buckling for an 
expansion of the cylinder. 


E2 


CRBL Circumferential stiffener buckling for a 

contraction of the cylinder. 

S.Y. Skin yield 

LRYU Longitudinal stiffener yield in tension 

LRYL Longitudinal stiffener yield in compression 

CRYU Circumferential stiffener yield in tension 

CRYL Circumferential stiffener yield in compres- 

sion 

WT Weight of the cylinder in pounds. 


The output obtained after each partial analysis is 
controlled by the user through the AESOP print control 
integers described later in this section. (See AESOP 
Namelist Input Description) . 


Namelist Data Block "CR1217" 


This data block is read and defined in the analysis 
subprogram NL1217. All nominal data values are established 
by the analysis subprogram D1217. Data block CR1217 
defines all the input variables required for the analysis 
subprograms. A complete list of the namelist data block 
CR1217 is presented in Table El. 


E3 


AESOP INPUT ROUTINES 


(o,o) 


m 


MAIN PROGRAM 

AESOP CONTROL PROGRAM AND PATTERN SEARCH 


ANALYSIS ROUTINES FOR SHELL PROGRAM 



FIGURE El. PROGRAM OVERLAY STRUCTURE 





TABLE El. 

NAMELIST Data Block "CR1217 


NMENOMIC 

NOMINAL 

VALUES 

DESCRIPTION 

BDV ( ) 

First seven 
true , eighth 
false. 

Logical variables determining 
active design variables. One 
for each design variable plus 
one to tell when the two depth 
variables are to be kept equal. 
The first seven quantities re- 
late one-for-one the seven 
design variables in the 
following order: t s t t X / t.<j) ,d x , 

d<j> , £x , £<|) * To make the two 
depth variables equal, the 
eighth value of the array BDV 
is made true, and the fifth 
value is made false . (The 

fifth value corresponds to d(j>) . 
The program will then make d<j> 
d x . 

BOTULX 

20.0 

i 

Upper limit for (d/t) x 

BOTULY 

20.0 

Upper limit for (d/t)^ 

CRCL ( ) 

-5. 0E4 

Compressive yield stress for 
circumferential stiffeners for 
each load condition, o, c/ ~„ 
<lbs/in2) <t,S0C 

1 

CRCU ( ) 

5.0E4 

, Tensile yield stress for circum- 
ferential stiffeners for each 
load condition , a^ S0T (lbs/in 2 ) 

DLT (2) 

0.0 

Indicator, zero when the longi- 
tudinal stiffeners are contin- 
uous; one otherwise, 6 . 

xw 

DLT (3) 

1.0 

Indicator, zero when circumfer- 
ential stiffeners are continu- 
ous,; one otherwise, 6^ . 


E5 



MNEMONIC 

NOMINAL 

VALUES 

DESCRIPTION 

EOF 

False 

Logical variable. If true, the 
program will terminate. 

EX( ) 

1.0E7 

Longitudinal stiffener modulus 
for each load condition, E xs . 

EY( ) 

1.0E7 

Modulus of circumferential 
stiffener for each load 
condition, 

El( ) 

1.0E7 

Longitudinal modulus of skin 
for each load condition E x 
(lbs/inch 2 ) 

E2( ) 

1.0E7 

Circumferential modulus of 
skin for each load condition, 

E.. 

4> 

Densities of the skin, circum- 
ferential stiffeners, and 
longitudinal stiffeners, re- 
spectively (lbs/inch2 ) 

GAM ( ) 

i .101 

GSM ( ) 

3750937.7 

Shear modulus of skin for each 
load condition, G. 

I 

1 

Number of load conditions, 
eight maximum; integer. 

ICACYC ( ) 

1, 2,3,4, 
5,6, If 8, 

9,10 

An array of ten elements used 
to specify optimization cycles 
on which to perform a complete 
analysis . 

I READ 

0 

Integer variable. If 1217, 
the program will read data 
cards in the same format as 
specified by NASA CR-1217. 


E6 






MNEMONIC 

NOMINAL 

VALUES 

DESCRIPTION 

IWRITE 

1217 

Integer variable. If 1217, 
the program will print the 
input data in the same format 
as shown in NASA CR-1217. 

KAPA ( , ) 

1.0 

Constants defining yield envel- 
ope, KaS- <TT is read first 
for each load condition then 
kct for each load condition. 
Similarly, kqq and ktc are read. 

L 

10.0 

Length of cylinder (inches). 

LRCL ( ) 

-5. 0E4 

Compressive yield stress for 
longitudinal stiffeners for 
each load condition, cr 
<lbs/in2) xS0C 

LRCU ( ) 

5.QE4 

Tensile yield stress for the 
longitudinal stiffeners for 
each load condition, a orvT1 
Ubs/in2) xS0T 

ML ( , ) 

20 

i 

Limit on the number of half 
wave numbers searched in the 
longitudinal direction for 
each load condition for each 
cylinder failure mode. The 
order is load condition then 
failure modes. Integer. 
(3x1). 

NL( , ) 

15 

Limits on the number of full 
wave numbers searched. 

NSM ( , ) 

(Not set) 

Number of modes saved for the 
approximate analysis. The 
values for the first load 
condition are read in the order 
gross, panel, skin, and then 
this is repeated for load 
condition two, etc. Integers 
(I x 3). 


E7 




MNEMONIC 

NOMINAL 

VALUES 

DESCRIPTION 

NUX ( ) 

.333 

Poisson's ratio of skin for 
each load condition, y x . 

NUY ( ) 

.333 

Poisson's ratio of skin for 
each load condition, 

NUl ( ) 

. 333 

Poisson's ratio for the circum- 
ferential stiffeners for each 
load condition- 

NU2 ( ) 

. 333 

Poisson's ratio for the longi- 
tudinal stiffeners for each 
load condition- 

Pl( ) 

0.0 

Applied axial compressive loads 
for each load condition, N 
(lbs/inch) . 

P2( ) 

0.0 

Applied external radial pressure 
for each load condition, p 
(lbs/inch 2 ) 

R 

10.0 

Radius of cylinder (inches) 

so( , ) 

l 

5.0E4 

Skin yield stresses. First 

S xOT rea< ^ f° r eac ^ load 

condition and then S x oc is 
read for each load condition. 
Similarly S<f>oT and S^oc are road. 


NOTE: Initial values of the design variables (t s ,tx,t<j) ,d x , 

d<j, , SL X > ) are input to the program through the AESOP namelist 

"IAESOP" as the ALPHA vector. 


E8 






NAMELIST Data Block "IAESOP 


This data block is read and defined in the optimization 
subprogram by subroutine BAESOP. All nominal data values are 
established by the optimization subroutine BDATA7. Data 
block IAESOP primarily defines which combination of the nine 
optimization search algorithms of the optimization subprogram 
are to be employed, how they are to be employed, and how many 
times the optimization cycle is to be repeated. 


The nine search algorithms available are described in 
reference 1; these are listed below. 


1. Sectioning 

2. Pattern 

3. Magnification 

4 . Steepest-Descent 


6 . Quadratic 

7 . Davidon 

8 . Random Point 

9 . Random Ray 


5. Adaptive Creeping 


A complete list of optimization data is presented in 
Tables E2. and E3 . . Table E2 contains the basic optimization 
control data. Table E3 contains the specialized print control 
data. It should be emphasized that all items in Tables E2 
and E3 are read by the single NAMELIST input block IAESOP. 


AESOP Print Control 

AESOP has a flexible print output capability. Varying 
levels of printout are available at user option as follows: 

Summary of function and control parameter values 
at the beginning and end of the optimization 
process ; 

Summary of function and control parameters values 
at the end of each cycle; 

Summary of function and control parameter values 
at the end of each evaluation; 

Detailed printout of individual search parameters. 

The convention adapted for print indicator is a six-letter 
mnemonic as described below. 

NOTE: In all cases print is given when the indicator is non- 

zero and omitted when it is zero. 
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EXAMPLES: 1) IPACRP=1, Print control vector following 

creeping search 

2) IPDPEN=1, Supply detailed print output 
from subroutine PENALTY 

IP X XXX Subroutine Reference 


CRP 

— 

CREEPR 

CYC 

= 

MAINOP (End of Cycle) 

DVD 

- 

DAVIDN 

MAG 

- 

MAGIFY 

PAT 

- 

PATERN 

PEN 

- 

PENLTY 

QUA 

- 

QUADRA 

RPT 

- 

RPOINT 

RRS 

- 

RANRAY 

SEC 

- 

SECTON 

STD 

- 

STDESC 


Type of Print 

A - Control Vector - ALPHA 

1 D - Detailed Print 

F - Function Array - FUNCTN 

' 1 Means PRINT INDICATOR 

An alphabetical list of print control indicators follows in Table E3. 
relevant search is identified for each input in the same way 
described earlier for the optimization data. 
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TABLE E2.. BASIC OPTIMIZATION DATA 


RELEVANT 

SEARCH 


MNEMONIC 123456789 


NOMINAL 

VALUES 


DESCRIPTION 


X 100*1* Determines first pertur- 

bation directions in 
creeping search. 

XXXXXXXXX 100*1. Nominal values of control 

parameters. 

XXXXXXXXX 100*1. Upper control parameter 

search limits. 

XXXXXXXXX' 100*1. Lower control parameter 

search limits. 


X XX 


X XX 


100 * 

.0000001 Minimum perturbations to 
be employed in creeping 
search. 

100*. 001 Starting perturbations 
for creeping search. 

.001 Initial termination tol- 
erance on Golden Section. 

.00001 Final termination toler- 
ance on Golden Section 


X X XX X XX X X 20*1. 


Final desired constraint 
tolerances . 


Steepest-descent weighting 
matrix indicator. 

0 - Unit matrix 

1 - Empirical matrix 

2 - Alternates between 

unit and empirical 
matrix 


Selects the order in 
which the control vari- 
ables are perturbed and 
sectioning searches. 

0 - Uniformly random 

1 - Natural order 

2 - Reverse natural order 


Used only in constraint logic 




RELEVANT 

SEARCH 


MNEMONIC 


IREPET 

ISECOF 



LIMIT X 


MAXCRP 


MAXDVD 


MAXJJJ 


X X 



MAXRPT 


MAXRRS 





NOMINAL 

VALUES 

DESCRIPTION 

10 

Number of optimization 
cycles to be completed. 

1000 

Cycle number at which 
sectioning search is 
terminated, . 

0 

Selects extreme of the 
search interval to be 
used when performance is 
constant on search ray. 


0 - Lower limit 

1 - Upper limit 

0 

Controls multiple extre- 
mal option. 


0 - Performance response 

surface unaltered 

1 - Performance response 

surface is warped 

2 

Number of sectioning 
searches. 

5 

Number of complete 
creeping searches to 
be performed when search 
is called. 

10 

The number of Davidon 
searches carried out when 
this method is selected. 

200 

The maximum number of 
performance evaluations. 

99 

Maximum number of magni- 
fication searches in an 
optimization calculation. 

10 

Number of random point 
evaluations . 

j 100 

i 

Number of random rays 
to be employed. 
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MNEMONIC 


RELEVANT 

S EARCH NOMINAL 

56(789 VALUES 


DESCRIPTION 


METHOP i XXXXXXXXX 1,2,3, 4, 5 


NALPHA X X XX X XX X X 3 


NFUNC XXXXXXXXX 1 


NMAXLO X X X X 10 


NMAXUP XX X X X 20 


NPHIAC XXXXXXXXX 1 

NPSIi* XXXXXXXXX 0 

NPTSRYt 10 

NUMOPT XXXXXXXXX 5 

NUMPSI* XXXXXXXXX 0 

NUMSTD X 2 

PHIEPS X XXX 0.0 


PSIWT . * XXXXXXXXX 20*. 

.0001 


The sequence of searches 
to be employed. 

1 - Sectioning 

2 - Pattern 

3 - Magnify 

4 ** Steepest-Descent 

5 - Creeping 

6 - Quadratic 

7 - Davidon 

8 - Random Point 

9 - Random Ray 

11- Arbitrary Ray Search 

Number of control param- 
eters to be employed 

Number of functions to be 
considered. 

Initial maximum number of 
evaluations in Golden 
Section . 

Final maximum number of 
evaluations in Golden 
Section and to limit the 
number of evaluations in 
pattern. 

Function number of the 
performance criteria. 

Constraint function numbers 

Number of steps to take 
when the ray search is 
called 

Number of optimization 
searches to be employed 

Number of constraints 

Number of steepest-desc- 
; ent searches 

! 

, Performance values with- 
in PHIEPS of the minimum 
value yet attained are 
treated as being equal 
in Golden Section 

Initial constraint error 
weights 


* Used only in constraint logic 
t Pertains to ray search 


El 3 






RELEVANT 

SEARCH 

NOMINAL 


MNEMONIC 

12 3 

4 5 6 

7 8 9 

VALUES 

DESCRIPTION 

QFACTR 


X 


1.0 

Quadratic perturbation 
factor 

QPERT^ 


X 


20*. 005 

Initial control param- 
eter perturbations for 
quadratic and Davidon 
searches 

RALOHIt 




True 

Logical variable used to 
indicate direction to go 
along the multidimensional 
ray. I.e., if true, go 
from XTENLO to XTENHI; if 
false, go from XTENHI to 
XTENLO. 

RANGEN 



X X 

1.0 

Random number generator 
trigger (1.0 means uni- 
form distribution) , (-1.0 
means normal distribution) 

RAYDIV+ 




10. 

Used to compute stepsize; 
i.e., if RAYDIV=10 , the 
stepsize will be such as 
to reguire 10 steps to go 
from XTENLO to XTENHI. 

RUFHI 



X 

' .01 

Maximum nondimens ional 
random ray perturbation 
size on any component. 
Value of 1.0 gives maxi- 
mum perturbation equal 
to search range. 

RUFLO 



X 

.00001 

Minimum nondimensional 
random ray perturbation 
size on any component 

SIBAR ± * 

XXX 

XXX 

XXX 

20*0.0 

Desired constraint values 

TOLFACj.* 

XXX 

XXX 

XXX 

20*0.5 

Constraint tolerance 
reduction factor 

TTOL^* 

XXX 

XXX 

XXX 

20*100. 

Initial constraint tol- 
erances 

WITER i * 

XXX 

XXX 

XXX 

100*1.0 

Starting values for iter- 
ative component of 
steepest-descent weighting 
matrix 


* Used only in constraint logic 
t Pertains to ray search 
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RELEVANT 

SEARCH 

NOMINAL 


MNEMONIC 

IXSl 

wm 

7 8 9 

VALUES 

DESCRIPTION 

WTDOWNi* 

X 

X X 

X 

20*0.5 

Constraint weight decrease 
factor 

XTENHI j_ 
XTENLOj_ 

XXX 

X 

XXX 
X X 

XXX 

X 

100*0.0 

100*0.0 

Used to extend upper search 
limit 

Used to extend lower search 
limit in Golden Section if 
performance is constant in 
feasible region 

i 

WARPALj_ 

XXX 

XXX 

IX X X 

100*0.0 

The warping origin in the 
control parameter space for 
multiple extremal feature 

WARPN 

XXX 

XXX 

1 

XXX 

2.0 

The exponent of the warping 
transformation 


*Used only on constraint logic 


El 5 














TABLE E3 . OPTIMIZATION PRINT CONTROL DATA 


RELEVANT 

SEARCH 


MNEMONIC t 1 2 3 


IPACRP 


IPACYC XXXXXXXXX 0 


IPADVD 


IPAMAG 


IPAPAT 


IPAPEN XXXXXXXXX 0 


IP AQUA 


IPARPT 


IPARRS 


IPASEC 


IPASTD 


IPDCRP 


IPDDVD 


NOMINAL 

VALUES 

DESCRIPTION 

0 

Creeping control param- 
eter print indicator. 

0 

Cycle control parameter 
print indicator. 

0 

Davidon control param- 
eter print indicator. 

0 

Magnification control 
parameter print indi- 
cator. 

0 

Pattern control param- 
eter print indicator. 

0 

Constraint penalty 
control parameter print 
indicator . 

0 

Quadratic control param- 
eter print indicator. 

0 

Random point control 
parameter print indi- 
cator. 

0 

Random ray control 
parameter print indi- 
cator. 

0 

Sectioning control 
1 parameter print indi- 
cator. 

0 

Steepest-descent control 
parameter print indi- 
cator. 

0 

Detailed creeping print 
indicator. 

0 

Detailed Davidon print 
indicator. 












MNEMONIC 


12 3 


RELEVANT 

SEARCH 

"IT'S” 6 |7 


IPDMAG 


X 


T 


8 9 


NOMINAL 

VALUES 


DESCRIPTION 


0 


Detailed magnification 
print indicator. 


IPDPAT 

IPDQUA 

IPDRIV 

IPDRPT 


X X 


X 


Detailed pattern print 
indicator . 

Detailed quadratic 
print indicator. 

Detailed derivatives 
print indicator. 

Detailed random point 
print indicator. 


IPDRRS 


IPDSEC 


IPDSTD 


IPFCRP 


X 


Detailed random ray 
print indicator. 

Detailed sectioning 
print indicator. 

Detailed steepest-de scent 
print indicator. 

Creeping function print 
indicator. 


IPFCYC 


XXX 


XXX 


XXX 


0 


IPFDVD 


X 


0 


Optimal function print 
indicator at the end of 
each cycle. 

Davidon function print 
indicator. 


IPFMAG X 

IPFPAT X 


0 

0 


Magnification function 
print indicator. 

Pattern function print 
indicator. 


IPFQUA 


X 


0 


Quadratic function print 
indicator. 


IPFRPT 


X 


0 


Random point function 
print indicator. 










MNEMONIC 


IPFRRS 


IPFSEC X 


IPFSTD 


IPGAIN 


IPNAML 


RELEVANT 

SEARCH I NOMINAL 

VALUES 



DESCRIPTION 


Random ray function print 
indicator. 

Sectioning function print 
indicator 

Steepest-descent function 
print indicator 

Print every iteration 
which will improve the 
performance. 

Namelist output control 
= 0, omit print 
= 1, print namelist data 
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AESOP Data Listings 


Program AESOP data can be conveniently grouped according 
to search and function. The user employing a particular 
search can independently specify the characteristics of that 
search and may not be concerned with input relevant to the other 
searches. Hence, a data grouping by search and by function is 
presented below for user convenience. It should be noted that 
certain inputs are common to more than one search; where this 
occurs, the input is repetitively defined in each search. 

Search Selection and Control . — 

NUMOPT - The number of optimization techniques to be employed. 
Each individual search request in a sequence of 
requests adds to this input (e.g., the search se- 
quence 4, 2, 4, 2 requires NUMOPT = 4). Maximum number 
of searches employed must satisfy NUMOPT < 20. 

METHOP. - The search sequence by numeric identification. For 
x example, the input METHOP (1) = 1, 2,3,4, 5,6,7, 8, 9 
signifies the following search sequence: 

1 - Sectioning 

2 - Pattern 

3 - Magnification 

4 - Steepest-Descent 

5 - Adaptive Creeping 

6 - Quadratic 

7 - Davidon (Fletcher-Powell) 

8 - Random Point (Monte-Carlo) 

9 - Random Ray (random evolution) 

11 - Arbitrary Ray 

The complete search sequence will be referred to as 
an optimization cycle. 

MAXJJJ - The maximum number of system evaluations. A direct 
iteration number limit. 

IREPET - The maximum number of times the search sequence (op- 
timization cycle) defined in METHOP^ will be utilized 

Parameter Selection . — 

NALPHA - The number of parameters available for optimization. 

No more than one hundred parameters may be employed. 
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ALPLCk - Lower bounds on each parameter search range 

ALPHI^ - Upper bounds on each parameter search range 

ALPHA. - The nominal parameter values. Note that ALPLO. < 

1 ALPHA* < ALPHI . must be satisfied. If a parti 1 

cular 1 parametei , / say ALPHA j , is to be fixed in 
value in a particular computation , then set 
ALPLOj = ALPHAj = ALPHI j . This effectively re- 
duces the parameter space dimension by one for 
each such parameter. 

Multiple Extremal Option . — 

IWARP - Controls multiple extremal option 

I WARP = 1, automatically warp the response surface 
IWARP = 0, leaves the response surface unmodified 

WARPAL. - The point at which the warping transformation is 
1 centered, i.e., the location of a known extremal 
point. 

WARPN - The degree of the warping transformation. The 
greater WARPN, the greater the response surface 
distortion. 

Optimization Function Selection . — 


FUNCTN. - AN INTERNAL ARRAY CONTAINING ALL COMPUTED OPTIMI- f 
1 Z AT ION FUNCTIONS j 


NFUNC - The total number of functions (FUNCTN ) being 
computed in the system model 
NOTE : NFUNC < 100. 

NPHIAC - The function to be minimized. AESOP always searches 
for a minimum; to maximize FUNCTN^ define FUNCTN n = 
-FUNCTN m and minimize FUNCTN n . 

NUMPSI - The total number of functions being constrained. 
NOTE: NUMPSI ^ 20. 

NPSI . - The functions to be constrained, e.g., NPSI(l) = 

3, 5, 1, 7 indicates that FUNCTN 3 , FUNCTN 5 , 

FUNCTNi, and FUNCTN 7 are to be constrained. 

SIBAR. - The desired values of the constraint functions 
1 defined by NPSH^ 
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FTOL ± 

TTOL i 

PSIWT ± 


WTUP. 

i 


WTDOWN. 

1 


The acceptable tolerances on the constraint 

function values, SIBAR. 

x 

Initial acceptable tolerances on the constraint 
function values, (should be approximately 100 
times greater than the corresponding FTOL^) . 

Initial constraint error weighting factors in the 
augmented performance function, <{>*, where 

4>* = <f> + (23) 

i 

Here 


W ± = PSIWT ± 

Incremental multiplicative constants used to 
increase the Wj_ on constraints which prove diffi- 
cult to satisfy. The nominal values of WTUPi = 
2.0 should be acceptable; hence, this input can 
normally be omitted. 

Decremental multiplicative constants used to 
decrease the Wj_ when a constraint is easily 
satisfied. The nominal values of WTDOWN^ = 0.5 
should be acceptable; hence, this input can 
normally be omitted. 


Sectioning Search Data (METHQP-j = 1) . — 

LIMIT - The number of times each parameter will be sec- 
tioned during a single sectioning search 

NMAXLO Maximum number of point evaluations employed in 

a single parameter's sectioning at search com- 
mencement (first optimization cycle) . In 
successive cycles, the maximum number of points 
employed is increased by one. 

NMAXUP An upper bound on the maximum number of point 

evaluations employed in sectioning a particular 
parameter . 

ISIDE - Indicator specifying selection of left or right 
boundary for a parameter that does not appear to 
affect the system performance 

ISIDE = 0, Select lower limits 
ISIDE = 1, Select upper limits 
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XTENHI. - Extension of higher search limits (ALPHIi) for a 
J ‘ parameter that does not appear to affect perfor- 
mance 

XTENLOj_ - Extension of lower search limits for a parameter 
that does not appear to affect performance 

IRANDM - Controls the order in which the parameters are 
sectioned 

IRANDM = 0, Random order selected 
IRANDM = 1, Natural order selected 
IRANDM = 2, Reverse natural order selected 

FACTHI - Section termination criteria. If three successive 
performance function values are within FACTHI 
of each other during sectioning of a given param- 
eter on the first optimization cycle, the section 
search of that parameter will cease. The termina- 
tion criteria is internally reduced with each 
optimization cycle. 

FACTLO - The lower limit on the termination criteria in 
any optimization cycle. 

ITRADE - Optimization/trade study indicator 

ITRADE = 0, Carry out a normal optimization search 
ITRADE = 1, Determine performance function sensi- 
tivity to each parameter by sectioning 
each parameter in turn about a given 
fixed point in parameter space. 

IPASEC, - Print indicators (See Table e 3, page E17) for section 

IPDSEC , search 

IPFSEC 

Pattern Search Data (METHQP-j = 2) . — 

IPAPAT, - Pattern Search print indicators 

IPDPAT , 

IPFPAT 


unification Search Data (METHOP-; = 3) . — 


MAXMAG - Maximum number of point evaluations performed 
during a single magnification search 
DELMAG - The magnification perturbation size, nominally 
set to 1% of distance to origin. Not normally 
modified from nominal value 
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IPAMAG, - Magnification search print indicators, (See Table 

IPDMAG, E3, Pages E16 to E18 ) . 

IPFMAG 


Steepest-Descent Search Data (METHOP^ = 4) . — 

NUMSTD - Number of gradient evaluations and one-dimensional 

searches performed each time that a steepest-descent 
search is requested during the optimization cycle. 

INDWMA - Steepest-descent weighting matrix indicator 
INDWMA = 0. , Unit matrix 
INDWMA = 1., Empirical matrix 
INDWMA = 2 • , Alternate on each cycle between 
unit and empirical matrices 

WITER^ - Learning factors for steepest-descent weighting 
matrix 


NMAXLO 


NMAXUP 


FACTHI 


FACTLO 


I PA STD , 
IPDSTD, 
IPFSTD 


Maximum number of point evaluations employed in 
the steepest-descent one-dimensional ray search at 
search commencement (first optimization cycle) . 

In successive cycles, the maximum number of point 
evaluations permitted is increased by one. 

Upper bound on the number of point evaluations 
along a steepest-descent one-dimensional ray in 
any optimization cycle. 

One-dimensional steepest-descent ray search 
termination criteria during first cycle. The 
termination criteria is reduced in each successive 
optimization cycle . 

Lower limit on one-dimensional steepest-descent 
ray search termination criteria, in any optimi- 
zation cycle. 

Steepest-descent search print control indicators. 


Adaptive Creeping Search Data (METHOPj = 5) - 

MAXCRP - Number of creeping search perturbations introduced 
into each parameter by a single adaptive creeping 
search in the optimization cycle. 
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IRANDM 


DCREEP i 

CREPMN i 

CREPMXi 

ALFSIN i 


I PAC RP , 
IPDCRP , 
IPFCRP 


Controls the order in which parameters are perturbed 

IRANDM = 0 , Random order 

IRANDM = 1 , Natural order 

IRANDM = 2 , Reverse natural order 

The initial perturbations to each parameter 

Minimum perturbations for each parameter 

Maximum perturbations for each parameter 

Direction of perturbation for each parameter 
(ALFSINj = ±1.0) 

Adaptive creeping search print indicators 


Quadratic Search Data (METHOPj = 6 ) 


QPERT ± 


Parameter perturbation magnitudes employed in 
computation of numerical partial derivative 
matrices r, 2 ^ 


3 oil 3 


and 


3<j> 

3Cd-j 


QFACTR 

NMAXLO 


NMAXUP 


FACTHI 


FACTLO 


IPAQUA, 
IPDQUA r 
IPFQUA 


- Scaling factor on the QPERT^ 

- Maximum number of point evaluations employed in 
the quadratic one-dimensional ray search at search 
commencement (first optimization cycle) . In suc- 
cessive optimization cycles , the number of point 
evaluations permitted increases by one. 

- Upper bound on the number of point evaluations 
along a quadratic one-dimensional ray search, in 
any cycle. 

- One-dimensional quadratic ray search termination 
criteria during first cycle. The termination 
criteria is decreased in each successive optimi- 
zation cycle. 

- Lower limit on one-dimensional quadratic ray 
search termination criteria, in any optimization 
cycle . 

- Quadratic search print indicators 
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Davidon Search Data (MET HOP j - 7) . 


MAXDVD 


QPERT ± 


NMAXLO 


NMAXUP 


FACTHI 


FACTLO 


IPADVD, 
IPDDVD , 
IPFDVD 


Number of Davidon (Fletcher-Powell) gradient 
evaluations and one-dimensional searches performed 
each time that a Davidon search is requested in 
the optimization cycle 

Parameter perturbation magnitudes employed in com- 
putation of numerical partial derivatives, 

3 <[> 
da j 

Maximum number of point evaluations employed in 
the Davidon one-dimensional ray search at search 
commencement (first optimization cycle) . In 
successive optimization cycles, the number of 
point evaluations permitted is increased by one. 

Upper bound on the number of point evaluations 
along a Davidon search one-dimensional ray in 
any cycle 

One-dimensional Davidon ray search termination 
criteria during first optimization cycle. The 
termination criteria is decreased in each succes- 
sive optimization cycle. 

Lower limit on one-dimensional Davidon ray search 
termination criteria, in any optimization cycle 

Davidon print control indicators* 


Random Point Search (MET HOP j = 8) . — 

MAXRPT - The maximum number of random points to be employed 
in the first request for a random point search 
within the optimization cycle. In successive 
requests, MAXRPT is set to zero, and no evaluations 
result. 

IPARPT, - Random point search print control indicators 
IPDRPT, 

IPFRPT 
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Random Ray Search (METHOPj = 9 ) 


MAXRRS 


RUFHI 


RUFLO 


IPARRS , 
PIDRRS , 
IPFRRS 


- The maximum number of random rays, one or 
two-sided, investigated each time the 
optimization cycle requests a random ray 
search. 

- The initial maximum non-dimension pertur- 
bation measure for each parameter. This 
is reduced each time random ray search 
consistently fails to improve performance. 

- Minimum-maximum dimensional perturbation 
measure for each parameter. 

- Random ray search print control indicators 


Arbitrary Ray Search (METHOP j = 11) . 

The arbitrary ray search searches the ray passing 
through two specified points in the multidimensional 
control space. For printout it is suggested that the 
tabular summary feature of AESOP be used. 


RALOHI 


NPTSRY 

RAYDIV 


XTENHIi 


Defines which direction to search? i.e., 
if TRUE, search "LO" + "HI", if 
FALSE, search "HI" + "LO" 

Total number of evaluations to be used 
on ray search 

Defines stepsize for ray search; i.e., 
control parameter step size = (XTENHIj_ - 
XTENLOi) /RAYDIV 

Defines the control parameter values 
at "HI" t end of the multidimensional 
ray to be searched. 


XTENHIi need not be greater than XTENLO^ , and XTENLOj 
need not be less than XTENHIj in this search . These two 
arrays merely define the end points of a ray in the multi- 
dimensional control space. 
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